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FINAL  TECHNICAL  REPORT 


Introduction 

The  principal  goal  of  this  research  effort  was  to  develop 
an  effective  computational  algorithm  for  the  estimation  of 
parameters  in  distributed  parameter  systems.  The  algorithm 
was  developed  in  a  general  setting  which  allowed  application 
to  phenomena  modeled  by  delay-differential  equations,  Volterra 
integral  equations,  and  partial  differential  equations  with 
memory  terms. 

In  particular  we  investigated  a  gradient-based  parameter 
estimation  method  for  dynamical  systems  in  an  abstract  space. 
The  basic  functional  analytic  framework  was  the  theory  of 
semigroups  of  operators  in  infinite  dimensional  space.  This 
framework  allowed  application  to  distributed  parameter  systems 
modeled  by  hereditary  systems  and  partial  differential 
equations.  This  research  focused  both  on  theoretical  aspects, 
such  as  convergence  criteria,  and  on  the  efficient 
implementation  and  testing  of  the  algorithms  for  computational 
purposes. 

Summary  of  research 

The  dynamical  systems  under  consideration  were  of  the 
general  form 


x(t)  =  4(q)xit)  +  Bu(t) 

x(0)  =  Xq  (1) 

y(t)  =  Cx(t) 

where  u  and  y  are  input  and  output  functions,  x  is  an  infinite 
dimensional  state,  and  j<(q)  is  an  evolution  operator 
depending  on  a  possibly  distributed  parameter  q.  A 
computationally  feasible  algorithm  was  developed  for  solving 
the  following  identification  problem. 

Problem  (ID)  .  Given  an  input  function  u  and  observgitions 
y^  at  times  t^,  i  =  1,  ...,  ra,  find  a  system  parameter  q 

which  minimizes  the  quadratic  cost  function 

m 

J(q)  *  E  iicx(t  ;q)  -  y  (1^ 

i  =  1 

where  x(t;q)  is  the  state  at  time  t  of  the  dynamical  system 
with  parameter  q. 
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Hereditary  models  of  fluid-structure  interaction  often 
contain  unknown  parameters  which  need  to  be  identified.  Such 
models  may  be  cast  as  a  dynamical  system  in  an  infinite 
dimensional  space  as  in  equation  (1).  Our  goal  was  to 
identify  delays  in  these  models  as  well  as  other  system 
parameters.  These  results  were  developed  jointly  with  J.  A. 
Burns  and  E.  M.  Cliff  of  Virginia  Tech  University  and 
published  in  a  paper  entitled,  "Parameter  identification  for 
an  abstract  Cauchy  problem  by  quasilinearization."  For 
immediate  reference,  this  paper  is  included  in  the  appendix  of 
this  report. 

Similarly,  certain  distributed  parameter  models  of 
viscoelastic  structures  may  also  be  formulated  as  an  abstract 
dynamical  system.  The  models  of  interest  contain  Boltzmann 
damping  and,  in  particular,  fractional  derivative  damping. 

The  numerical  phase  of  this  research  effort  was  principally 
concerned  with  the  following  aspects  of  the  general  meeting 
described  above; 

(1)  the  development  and  testing  of  a  numerical  algorithm  for 
estimating  parameters  in  a  Volterra  integral  equation 
arising  from  a  viscoelastic  model  of  a  flexible  structure 
with  Boltzmann  damping; 

(2)  the  implementation  of  numerical  methods  for  a  system  of 
Volterra  equations  resulting  from  a  Galerkin 
approximation  of  a  partial  differential  equation  with 
hereditary  effects. 

Our  principal  research  effort  was  directed  toward  the 
development  and  estimation  of  distributed  parameter  models  of 
flexible  structures  with  internal  damping.  The  design  of 
control  systems  for  flexible  structures  is  highly  dependent  on 
the  amount  of  internal  damping  present  in  the  structure. 
Damping  parameters  typically  change  as  materials  and 
geometries  of  the  structures  change.  Accurate  and  efficient 
identification  algorithms  are  needed  to  estimate  a  system's 
characteristics  and  implement  a  stable  control  algorithm. 

In  particular,  we  considered  the  partial  differential 
equation 

pu^^(x,t)  =  |Eu^(x,t)  +  J^g{t-s)Uj^(x,s)dsJ. 

+  f(x,t)  (2) 

onO<x<l,  t>0,  with  appropriate  boundary  and  initial 
conditions.  The  function  u(x,t)  represents  the  longitudinal 
displacement  at  position  x  and  time  t  along  a  uniform  bar  of 
density  p  where  E  is  a  stiffness  parameter  and  f(x,t)  is  a 
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forcing  function.  The  function  g(s)  models  damping  effects. 
The  integral  represents  a  Boltzmann-type  internal  damping  term 
which  assumes  that  the  stress  is  function  of  the  strain  and 
the  strain  history.  We  studied  fractional  derivative  damping 
models  in  which  the  kernel  function  is  given  by  an  expression 
of  the  form 


g(s)  =  — 2 - -  ,  s  >  0,  0  >  0,  0  <  a  <  1. 

r(l-a)s“ 

An  estimation  algorithm  for  a  discretized  form  of  equation  (2) 
was  formulated  and  tested.  Using  simulated  data,  it  was  shown 
to  be  means  of  identifying  the  parameters  a  and  0. 

A  Galerkin  approximation  of  equation  (2)  using,  for 
example,  linear  splines  in  the  space  variable  yields  a  system 
of  Volterra  integro-dif ferential  equations  with  an  integrable, 
but  unbounded,  kernel.  A  gradient-based  identification 
algorithm  for  the  parameters  a  and  (3  was  implemented  in  a 
Volterra  equation  of  the  form 

’  w(t)  =  Mw(t)  +  j  K(-s)w(t+s)ds  +  F(t),  t  >  0, 

.  W(0)  =  TJ,  W(S)  =  ^(S),  S  <  0, 


This  involved  the  implementation  of  numerical  methods  for 
solving  the  Volterra  equation  and  its  sensitivity  equations 
with  respect  to  the  unknown  parameters.  Where  possible, 
numerical  results  were  checked  against  a  closed-form  solution 
obtained  by  a  Laplace  transform  method  using  software  capable 
of  symbolic  computation. 

A  gradient-based  algorithm  for  identifying  a  singularity 
in  a  weakly-singular  Volterra  integral  equation  was 
established  and  numerically  tested.  These  results  were 
published  in  the  Journal  of  Integral  Equations  and 
Applications  in  a  paper  referenced  below.  Additional 
numerical  results  were  published  in  Applied  Numerical 
Mathematics ,  also  referenced  below.  For  the  immediate 
reference  this  paper  is  included  in  the  appendix  together  with 
several  figures  which  did  not  appear  in  the  published  version. 

The  Galerkin  approximation  of  the  hereditary  partial 
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differential  equation  model  in  the  space  variable  gives  rise 
to  a  large  system  of  Volterra  equations  with  weakly  singular 
kernels.  Our  experience  indicated  the  importance  of  an 
accurate  and  efficient  method  for  solving  systems  of  this  type 
with  non-smooth  solutions.  We  employed  a  product  integration 
method  in  the  time  variable. 

A  class  of  fractional  linear  multi-step  methods  for  the 
numerical  solution  of  weakly  singulary  Volterra  equations  was 
investigated,  but  not  found  to  be  significantly  superior  to 
product  methods  in  this  context.  This  was  attributed  to  the 
fact  that  linear  multi-step  methods  require  a  priori  knowledge 
of  the  singularity.  Numerical  results  indicated  that  the 
convergence  rates  of  these  methods  deteriorate  in  the  absence 
of  this  knowledge.  In  our  context  the  singularity  is  among 
the  parameters  to  be  identified  so  is  not  known  a  priori.  Our 
problem  required  a  numerical  integration  scheme  which  is 
robust  over  a  wide  range  of  singularities. 

The  final  numerical  results  of  this  research  effort 
concerned  a  method  for  parameter  estimation  in  a  semi-discrete 
approximation  of  the  partial  differential  equation 


pu^^(x,t)  =  Eu^^(x,t) 


t  -0(t-s) 

f 


r{i-a)  at  '0  (t-s) 


with  boundary  conditions  u(0,t)  =  0,  u(l,t)  =  0,  and 
initial  conditions  u(x,0)  =  Uq{x),  u^(x,0)  =  Uj^(x)  . 
Integrating  with  respect  to  t  one  obtains 


r  pu^(x,t) 


=  e|  Uj^j^(x,s)ds  + 

I 


r(l-a)  •'0  (t-s) 


t  -/3(t-s) 

- 5r 


where  f^(x,t) 


r  f  (x,s)ds  +  pu  (X) . 

J  0 


We  applied  a  Galerkin  approximation  in  which  the  interval 
[0,1]  is  divided  into  N  equal  parts  and  the  homogeneous 
boundary  conditions  allowed  to  approximate  the  solution  by  a 

N 

function  of  the  form  u(x,t)  =  V  a.(t)(^,  (x)  where  is  a  cubic 

j=0  ^  J  ^ 

spline  basis  element.  Substituting  in  the  equation,  taking  an 
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inner  product  with  a  basis  element  (p^,  and  integrating  by 
parts  yielded  a  Volterra  equation  in  t  of  the  form 


pA^v(t) 


e“/3(t-s) 

(t-s)“ 


v(s) ds 


f3(t) 


N  N  .  . 

where  A  and  D  are  N+1  x  N+1  matrices  depending  on  inner 

products  of  0 j  and  <P'^ . 

The  quasilinearization  algorithm  required  that  we  solve 
this  equation  along  with  its  sensitivity  equations  obtained  by 
differentiating  with  respect  to  a  and  /3,  the  unknown 
parameters.  The  sensitivity  equations  also  have  weakly 
singular  kernels  so  that  the  same  numerical  methods  may  be 
applied.  These  results  will  be  reported  in  the  Proceedings  of 
the  World  Congress  of  Nonlinear  Analysts .  This  paper  in 
included  for  reference  in  the  appendix  of  this  report. 


Research  articles 

The  following  papers  relating  to  this  research  effort  were 
published  or  will  soon  appear  in  refereed  journals  or  as 
invited  papers  in  conference  proceedings. 

D.  W.  Brewer,  Gradient  methods  for  identification  of 
distributed  parameter  systems.  Proceedings  of  the  28th  IEEE 
Conference  on  Decision  and  Control,  December  1989,  599-603. 

D.  W.  Brewer  and  R.  K.  Powers,  Parameter  identification  in  a 
Volterra  equation  with  weakly  singular  kernel.  Journal  of 
Integral  Equations  and  Applications  2(1990),  353-373. 

D.  W.  Brewer  and  R.  K.  Powers,  Parameter  estimation  for  a 
Volterra  integro-dif ferential  equation.  Applied  Numerical 
Mathematics  9(1992),  307-320. 

D.  W.  Brewer  and  R.  K.  Powers,  A  direct  method  for  parameter 
estimation  in  distributed  systems.  Proceedings  of  the  World 
Congress  of  Nonlinear  Analysts ,  August  1992,  to  appear. 

D.  W.  Brewer,  J.  A.  Burns,  and  E.  M.  Cliff,  Parameter 
estimation  for  an  abstract  Cauchy  problem  by 

quasilinearization.  Quarterly  of  Applied  Mathematics  51(1993), 
1-22. 
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PARAMETER  IDENTIFICATION  FOR  AN  ABSTRACT  CAUCHY  PROBLEM 


BY  QUASILINEARIZATION 

Dennis  W.  Brewer 

University  of  Arkansas,  Fayetteville 

John  A.  Burns 

Virginia  Polytechnic  Institute  and  State  University 
Eugene  M.  Cliff 

Virginia  Polytechnic  Institute  and  State  University 

Abstract.  A  parameter  identification  problem  is  considered  in  the  context 
of  a  linear  abstract  Cauchy  problem  with  a  parameter-dependent  evolution 
operator.  Conditions  are  investigated  under  which  the  gradient  of  the 
state  with  respect  to  a  parameter  possesses  smoothness  properties  which 
lead  to  local  convergence  of  an  estimation  algorithm  based  on  quasi¬ 
linearization.  Numerical  results  are  presented  concerning  estimation  of 
unknown  parameters  in  delay-differential  equations. 


The  research  of  the  first  and  second  authors  was  supported  in  part  by  the 
National  Aeronautics  and  Space  Administration  under  contract  NASl-18107 
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Applications  in  Science  and  Engineering  (ICASE),  NASA  Langley  Research 
Center,  Hampton,  VA  23665.  The  research  of  the  first  author  was  supported 
in  part  by  the  Air  Force  Office  of  Scientific  Research  under  grant 
AFOSR-89-0472 .  Additional  support  was  received  from  DARPA  under  grant 
F49620-87-C-0116  while  the  first  author  was  in  residence  at  the 
Interdisciplinary  Center  for  Applied  Mathematics,  Blacksburg,  VA  24061. 
Also,  research  of  the  second  and  third  authors  was  supported  in  part  by 
APOSR  under  grant  APOSR-89-0001 ,  DARPA  under  grant  F49620-87-C-0116,  and 
SDIO  under  contract  F49620-87-C-0088. 
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1.  Introduction 


During  the  past  fifteen  years  considerable  effort  has  been  devoted  to 
the  problem  of  estimating  unknown  parameters  in  distributed  parameter 
systems.  The  recent  book  by  Banks  and  Kunisch  [9]  provides  an  excellent 
account  of  the  progress  made  in  the  field.  Many  parameter  estimation 
problems  are  best  formulated  as  optimization  problems  (sometimes  over 
infinite  dimensional  "parameter  spaces”)  and  algorithms  are  developed  to 
minimize  an  appropriate  cost  function.  Although  there  are  several 
approaches  to  these  problems,  their  infinite  dimensional  nature  requires 
that  numerical  approximations  be  introduced  at  some  point  in  the  analysis. 
Consequently,  there  are  two  basic  classes  of  algorithms  for  optimization 
based  parameter  estimation.  The  first  type  of  algorithm,  and  the  most 
frequently  used  for  dynamic  problems,  is  indirect  and  proceeds  by  initially 
approximating  the  dynamic  equations  (e.g.  finite  elements,  finite 
differences,  etc.)  and  then  using  optimization  algorithms  on  the  finite 
dimensional  problem.  This  approach  is  typified  by  the  papers  [1]~[6],  [8], 
[10],  and  [18]. 

The  second  more  direct  approach  is  based  on  the  direct  application  of 
an  optimization  algorithm  and  employing  numerical  approximations  at  each 
step  of  the  algorithm  to  compute  the  necessary  solutions  of  the  dynamic 
equations.  This  approach  is  used  in  [12],  [13],  [17],  and  [19].  Both 
methods  have  advantages  and  disadvantages.  Depending  on  the  particular 
type  of  distributeo  parameter  system,  one  method  may  out  perform  the 
other . 

Although  we  shall  consider  only  the  problem  of  identifying  a  finite 
number  of  parameters,  the  infinite  dimensional  dynamic  constraint  enters 
into  the  optimization  algorithm.  Basically,  the  objective  function  from 
parameter  space  to  R  is  a  composition  of  a  finite  rank  map  with  an 
operator  (defining  the  dynamic  constraint)  on  an  infinite  dimensional 
space.  Therefore,  any  method  that  requires  gradients  to  be  computed  will 
have  to  deal  with  the  differentiation  of  the  infinite  dimensional 
constraint,  i.e.  the  chain  rule  is  needed.  It  is  in  this  sense  that  the 
quasilinearization  algorithm  considered  here  has  an  "infinite  dimensional" 
nature . 
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Direct  methods  such  as  quasi  1 inearization  considered  here  are  often 
limited  by  the  fact  that  the  dependence  on  unknown  parameters  of  the 
solution  to  the  infinite  dimensional  dynamical  equations  may  not  be  "smooth 
enough"  to  establish  convergence  of  the  algorithm.  Indeed,  some  algorithms 
may  not  be  properly  defined  without  this  necessary  smoothness.  Indirect 
methods  avoid  this  difficulty  and  often  lead  to  easily  implemented 
algorithms.  On  the  other  hand,  when  direct  methods  can  be  applied  it  is 
sometimes  possible  to  establish  the  convergence  and  the  rates  of 
convergence  to  the  unknown  optimal  para.Ticters  (see  [13],  [19]). 

This  paper  considers  the  dependence  on  an  unknown  parameter  q  of  the 
solution  of  the  linear  abstract  Cauchy  problem 


(  x(t)  =  A(q)x(t)  +  u(t).  0  <  t  <  T, 

(x(0)=x„. 

Our  ultimate  g)al  is  to  formulate  and  establish  the  convergence  of  a 
gradient-based  parameter  estimation  algorithm  applicable  in  this  abstract 
setting. 

This  algorithm  employs  computation  of  the  gradient  D^x(t;q)  of  the 

solution  of  (1.1)  with  respect  to  the  parameter.  Conditions  for  the 
existence  of  this  gradient  are  established  in  [11],  In  Section  2  we  review 
these  conditions  and  the  general  setting  for  the  remainder  of  the  paper. 
Convergence  of  the  algorithm  requires  certain  smoothness  properties  of  the 
gradient  D^x(t;q)  with  respect  to  q.  These  properties  are  established  in 

Section  3  and  their  applicability  to  a  linear  delay-differential  equation 
is  discussed  in  Section  4.  In  this  example  the  delay  is  among  the 
parameters  so  that  in  this  setting  the  parameter  dependence  appears  in 
unbounded  terms  of  the  evolution  operator  A(q). 

An  abstract  parameter  estimation  algorithm  for  a  finite  dimensional 
parameter  space  using  a  discrete  cost  function  is  presented  in  Section  5. 

In  Section  6  its  convergence  is  established  using  the  results  of  Section  3. 
In  Section  7  we  present  several  numerical  examples  which  indicate  the 
performance  of  the  algorithm  for  delay  and  coefficient  estimation  in  linear 
delay-differential  equations.  Additional  examples  may  be  found  in  [I2l. 
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Numerical  testing  and  evaluation  on  a  wider  variety  of  parameter  estimation 
problems  will  be  undertaken  in  a  subsequent  paper. 

2.  The  General  Setting 

The  application  of  quasilinearization  to  parameter  estimation  requires 
knowledge  of  the  derivative  of  the  state  with  respect  to  the  unknown 
parameter.  This  topic  is  addressed  in  [11].  In  this  section  we  review  the 
framework  used  there  to  obtain  differentiability  and  establish  notation  to 
be  used  in  the  remainder  of  this  paper. 

Let  P  be  an  op^n  subset  of  a  normed  linear  space  P  with  norm  I*]  and 
let  X  be  a  Banach  space  with  norm  jj•jj.  For  every  q  e  P  let  A(q)  be  a 
linear  operator  on  D(A(q))  in  X.  Throughout  this  paper  we  assume 

(HI)  A(q)  generates  a  strongly  continuous  semigroup  S(t;q)  on  X; 

(H2)  D(A(q))  =  D  is  independent  of  q; 

(H3)  ||S(t;q)xl|  <  Me^^||x||,  x  €  X,  t  >  0,  q  €  D,  for  some  constants 

M  and  a/  independent  of  q,  x,  and  t. 

i 

Fix  T  >  0  and  u  e  L  (0,T;X).  Define  Q(t;q)  =  S( t-s ;q)u{ s)ds  for  q  6  P, 

0  <  t  <  T.  Note  that  if  (1.1)  has  a  strong  solution  then  it  is  given  by 

the  formula  x(t)  =  S(t;q)xQ  +  Q(t;q)  for  0  <  t  <  T. 

In  applications  of  this  theory  it  is  useful  to  consider  just  those 
terras  of  A(q)  in  which  the  parameter  appears.  To  this  end  we  write 
A(q)  =  A  +  B(q)  where  A  and  B(q)  both  have  domain  D  and  A  is  independent 
of  q.  Concerning  B(q)  we  assume  the  following; 

(H4)  For  every  q,  e  P  there  is  a  constant  K  such  that 

T 

r  ||B(q)S( t;q„)x||dt  <  Kllx||  for  all  x  €  D. 

Jo  ^ 
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In  Section  4  we  discuss  an  example  in  which  an  unbounded  operator  B(q) 
satisfies  (H4) .  This  hypothesis  does  imply,  however,  that  the  linear 

mapping  x  ->  B(q)S(*;qQ)x  is  bounded  as  a  mapping  from  D  into  L*(0,T;X). 

Let  Flq.q^)  denote  the  bounded  linear  extension  of  this  operator  to  X.  Let 

denote  the  norm  in  L^(0,T;X).  Concerning  F  we  assume  the  following; 

(H5)  There  is  closed  subspace  Y  of  X  such  that 

(i)  F(q,qQ)Xg  €  lNo,T;Y)  for  q,  q^  6  P,  and 

(ii)  for  every  €  P  and  f  >  0  there  exists  S  >  0  such  that 
llF(q,qQ)y  -  for  y  €  Y  and 

|q  -  qol  ^ 

The  analogue  of  F  for  the  function  Q(t;q)  is  the  mapping  GCq.q^)  from 
lVo,T;D)  into  L^(0,T;X)  defined  by 

tGlq.q^lwK t)  =  J  B{q)S(t-s;qQ)w(s)ds. 

By  (H4)  is  follows  that  G  can  be  extended  to  a  bounded  linear  mapping  on 
L^(0,T;X)  so  that  in  particular  GCq.q^lu  is  defined  as  an  element  of 

L*(0,T;X).  In  addition  we  assume 

(H6)  GCq.q^lu  €  L^0,T;Y)  for  q.  q^  €  P 

where  Y  denotes  the  subspace  required  by  (H5). 
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3.  Parameter  Dependence 


In  this  section  we  deduce  smoothness  properties  of  the  solution 
x(t:q)  =  S(t;q)xQ  +  Q(t;q)  with  respect  to  q.  These  properties  are  derived 

from  similar  properties  of  F(q,qQ)  and  G(q,qQ)  which  are  operators  related 

to  A(q).  These  results  will  be  used  in  Section  5  to  prove  convergence  of 
the  parameter  estimation  algorithm.  Throughout  this  section  T  >  0,  x^  e  X, 

and  u  €  L^(0,T;X)  are  fixed  as  given  in  (1.1).  The  symbol  denotes 

Frechet  differentiation  with  respect  to  q.  These  results  are  given  as  a 
series  of  lemmas  whose  proofs  are  at  the  end  of  this  section. 


Lemma  3,1.  Suppose  (HI)  -  (H5)  hold.  In  addition,  suppose  that  for  a 
given  q*  €  P 


(H7)  F(q,qQ)xQ  is  Frechet  differentiable  with  respect  to  q  at  q^ 
for  every  q^  €  P. 


For  brevity,  let  DF(qQ)  denote  D  tf’(q,qQ)xQ]  |  for  q^  e  P.  In  addition 


suppose 


(H8)  DF(q)  is  strongly  continuous  in  q  at  q*,  that  is,  for  each 

h  e  P  the  mapping  q  DF(q)h  from  P  into  L^(0,T;X)  is 
continuous  at  q*. 


Then  for  each  t  e  [0,T],  S(t;q)xQ  is  Frechet  diffent  table  with  respect  to  q 
at  every  q  e  P  and  D^[S(t;q)xQ]  is  strongly  continuous  with  respect  to  q 
at  q* . 

Lemma  3.2.  Suppose  (HI)  -  (H6)  hold  and  in  addition  suppose  that  for  a 
gi ven  q*  e  P, 
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(H9)  G(q,qQ)u  is  Frechet  differentiable  with  respect  to  q  at  q^ 
for  every  q^  €  P. 

Again  denoting  this  derivative  by  DCCq^)  for  q^  6  P,  assume 
(HIO)  DG(q)  is  strongly  continuous  in  ct  at  q*. 

Then  for  t  €  [0,T],  Q(t;q)  is  Frechet  differentiable  with  respect  to  q  at 
every  q  e  P  and  D^[Q(t;q)]  is  strongly  continuous  in  q  at  q* . 

Lemma  3.3.  Suppose  (HI)  -  (H5)  and  (H7)  hold  and  in  addition  suppose 

(Hll)  F(q,q*)  is  locally  Lipschitz  continuous  in  q  at  q* ,  uniformly 

for  y  €  Y,  that  is,  there  exist  constants  >0  such  that 

||F(q,q*)y  -  Frq.q*)y||^  <  Kjq  -  q*l  lly|| 
whenever  |q  -  q*|  <6^  and  y  €  Y. 

Moreover,  assume  that 

(H12)  DF(q)  is  strongly  locally  Lipschitz  continuous  with  respect 

to  q  ar  q*.  That  is,  for  each  h  €  P,  there  are  constants 
K,  6  >  0  such  that 

|iDF(q)h  -  DF(q*)hl|  <  "  q*! 

for  |q  -  q*l  <6. 

Then  D^[S(t;q)xQ]  is  strongly  locally  Lipschitz  continuous  with  respect  to 
q  at  q*  for  every  t  e  [0,T]. 
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Lemma  3.4.  Suppose  (HI)  -  (H6),  (H9)  -  (HIO)  hold  and  in  addition  suppose 


{H13)  DG(q)  is  strongly  locally  Lipschitz  continuous  with 
respect  to  q  at  q*. 

Then  Dg[Q(t;q)]  is  strongly  locally  Lipschitz  continuous  with  respect  to  q 
at  q*  for  every  t  e  [0,Tl. 

Although  the  assumptions  (Hi)  -  (H13)  are  rather  technical,  we  shall 
see  that  they  can  be  easily  verified  for  delay  systems  even  in  the  case 
that  the  unknown  parameter  is  the  delay  itself.  Therefore,  the  results 
presented  here  remove  the  limitations  placed  on  the  perturbation  B(q)  in 
papers  [13]  and  [161. 

For  completeness  we  now  present  the  proofs  of  Lemma  3.1  -  Lemma  3.4. 
However,  these  proofs  make  use  of  the  basic  results  found  in  [11]  and  in 
order  to  keep  the  length  of  the  proofs  reasonable  we  assume  that  the  reader 
has  [11]  in  hand. 

Proof  of  Lemma  3.1.  It  is  shown  in  [11]  that  (HI)  -  (H5),  (H7)  imply  that 
D^[S(t;q)xQ]  exists  for  q  6  P.  Furthermore,  it  is  given  by  the  formula 

(3.1)  D^[S(t;q)x_]h  =  f*S(t-s;q)[DF(q)h](s)ds,  h  €  P. 

q  ^  Jq 

We  therefore  obtain  by  substitution 

(3.2)  D^[S(t;q)xQ]h  -  Dq[S( t ;q*)xQ]h 

=  f  [S(t-s;q)  -  S(t-s;q*))([DF(q)h](s))ds 

+  f^S(t-s;q*)([DF(q)h](8)  -  [DF(q*)h] (s) )ds . 

"^0 
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Let  f  >  0  be  given  and  let  C  =  Me^^.  It  can  be  shown  (see  the  proof  of 
Theorem  1  [11])  that  for  all  x  e  X 

(3.3)  l|S(t;q)x  -  S(t:q*)x||  <  Cl|F(q.q*)x  -  F(q*,q*)xl|j . 

Combining  (3.3)  with  (H5ii)  shows  that  for  some  >  0 

||S(t,q)y  -  S(t;q*)y((  <  £C(|y||,  0  <  t  <  T,  y  e  Y, 

whenever  |q  -  q*|  <  ^ ^ •  In  particular,  putting  y  =  [DF(q)h](8)  €  Y  by 
(H5i)  we  obtain 

|l[S(t-s;q)  -  S(t-s;q*)][DF(q)h](s)l|  <  cCl|[DF(q)h]( s ) H 

for  |q  -  q*|  <  6^,  a.e.  s  6  (0,T).  Since  DF(q)h  is  continuous  at  q* ,  there 

exist  constants  K  ,  6,  >0  such  that 
2  2 

||DF(q)h||j  <  Kj  for  |q  -  q*|  <  ^ ^ . 

Combining  these  estimates  shows  that  the  first  term  in  (3.2)  is  bounded 
by  if  |q  -  q* |  <  min(6^,5j). 

Using  (H8)  it  is  easy  to  see  that  there  exists  6  >  0  such  that  the 

O 

second  term  in  (3.2)  is  bounded  by  eC  for  jq  -  q*l  <  6  .  These  estimates 

3 

complete  the  proof  of  Lemma  3.1. 
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Proof  of  Lemma  3.2.  By  Theorem  3  of  [11],  D^[Q(t;q)]  exists  for  q  e  P  and 

(3.4)  D^[Q(t;q)]  -  Dq[Q(t;q*)] 

=  f  [S(t-s;q)  -  S(t-3;q*)][DG(q)(s)]ds 

+  f^S(t-s;q*)[(DG(q))(s)  -  (DG(q*) ) ( s ) Ids 
‘'O 

where  u  has  been  suppressed  in  the  notation.  Since  IXj(q)  e  L^(0,T;Y)  for 
q  6  P  by  (H6),  the  proof  follows  exactly  as  in  the  proof  of  Lemma  3.1. 

Proof  of  Lemma  3.3.  Let  £  >  0  be  given.  By  (3.3)  and  (Hll)  there  exists 
>  0  such  that 

|lS(t;q)y  -  S(t;q*)yl|  <  CKjly|||q  -  q*| 

for  y  6  y  and  jq  -  q*|  ^  •  Since  DF(q)h  €  L^(0,T;Y)  by  (Hoi)  we  have  as 

in  the  proof  of  Lemma  2.1  that  the  first  term  of  (3.2)  is  bounded  by 

K  K  Iq  -  q*(  for  |q  -  q*|  <  min  (^,,^,).  An  estimate  of  the  same  form  is 

easily  obtained  for  the  second  term  of  (3.2)  using  (H12).  These  estimates 
complete  the  proof  of  Lemma  3.3. 

Proof  of  Lemma  3.4.  Since  DG(q)u  €  L^(0,T;Y)  by  (H6),  the  proof  follows 
exactly  as  in  the  proof  of  Lemma  3.3  using  (3.4)  in  place  of  (3.2). 

4.  Application  to  a  Delay-Differential  Equation 

In  this  section  we  apply  the  framework  of  the  previous  sections  to  the 
linear  delay-differential  equation 
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(4.1) 


x{t)=  af.x(t)  +  E  a,x(t 
^  k=l 
x(0)  =  r\ 


-  qj^)  +  u(t) 


Let  P  =  r”,  fix  r  >  0,  and  let  P  =  {q  =  (q^.q^,  .  .  .  ,  q^^)  :  0  <  qj^  <  r 
for  k  =  1,2,.  .  .,n}.  In  equation  (4.1),  fj  e  JR,  e  R,  k  =  0,1,.  .  .,n, 

tp  6  L^(-r,  0)  with  norm  denoted  by  jl^||  ,  u  e  L^(0,T),  and  x.is)  =  x(t+s) 

1  V 

for  t  >  0,  -r  <  s  <  0.  By  a  solution  of  (4.1)  we  mean  a  function  x  which 
is  absolutely  continuous  on  [0,T]  and  satisfies  (4.1)  almost  everywhere  on 
(0,T). 

Following  the  construction  in  [14],  we  take  X  =  R  x  lV-t.O)  with  norm 
l|('7. 9^)11  =  l»7l  +  and  define  for  q  €  P  an  operator  A(q)  on 


D  =  {(r},'p)  e  X:  ^  is  abs.  cont.  on  [-r,0],  p  €  L^(-r,0),  and 

»3(0)  =  Tj) 


11 

A(q)(T7,p)  =  (  ,  V  ) 


Then  is  well  known  that  A(q)  generates  a  strongly  continuous  semigroup 
S(t;q)  on  X  satisfying  S(t;q)  =  (y(t),  y^)  where  y(t)  =  y(t;q)  denotes  the 

solution  of  (4.1)  with  u  =  0.  It  is  a  consequence  of  standard  results  that 
(Hi)  -  (H3)  hold  in  this  setting. 

For  q  =  (q^,.  .  idjj)  and  q^  in  P,  (>?,»?)  e  X,  and  w  e  L^(0,T)  it 

follows  that  in  this  example  the  mappings  F  and  G  of  Section  3  are  given  by 
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(4.2) 


F(q,qQ)(»7,^)  =  ( 


n 

lc~l 


0  ) 


and 


(4.3) 


(G(q,qQ)w}( t)  =  ( 


n 


for  a.e.  t  e  (0,T)  where  z(t;q)  denotes  the  solution  of  (4.1)  with  u  =  w 
and  ir],ip)  =  (0,0).  It  is  shown  in  [11]  that  these  mappings  satisfy 
(H4)  -  (H6)  with  the  closed  subspace  Y  =  R  x  {0}.  It  is  also  shown  in  [11] 
that  F  and  G  satisfy  the  differentiability  hypotheses  (H7)  and  (H9)  for 
=  Xq  e  D  and  q.q^  e  P-  Furthermore,  their  Frechet  derivatives  are 

given  by 


(4.4) 


[DF(q)h](t)  =  (  -  E  aj^y(t-qj^;q)hj^ 
k*l 


0  ) 


and 

n 

(4.5)  [DG(q)h](t)  =  (  -  E  aj^z( t-qj^;q)hj^,  0  ) 

k*l 


for  q  €  P,  h  =  (hj . h^)  e  R*',  where  y(t;q)  is  the  solution  of  (4.1) 

with  u  =  0  and  z(t;q)  is  the  solution  of  (4.1)  with  (»7,y?)  =  (0,0). 

It  remains  to  establish  conditions  under  which  (H8),  (HlO)  -  (H13)  are 
satisfied. 

Lemma  4.1.  Fix  q*  =  (q*,.  .  •  .q*)  €  P  and  e  D.  Then  F(q,q*)xQ  as 
defined  by  (4.2)  satisfies  (Hll). 

Proof:  In  Section  5  of  [11]  it  is  shown  that  there  are  constants  C  and 
-  2 
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6  >  0  such  that 

2 


llF(q*+h,q*)(r7,0)  -  F(q* ,q*) jh| |l(rj,0)|| 


n  " 

for  h  €  R  ,  »7  e  R,  |h|  <  6_.  Here  we  define  |h|  =  E  |h,  |  .  This  estimate 


k=l 


is  equivalent  to  (Hll)  with  Y  =  R  x  {0}. 


Lemma  4.2.  Suppose  e  D.  Then  DF(q)  as  given  by  (4.4)  satisfies 

(H8) .  Moreover,  if  in  addition  tp  is  of  bounded  variation  on  [-r,0],  then 
DF(q)  satisfies  (H12) . 

Proof :  Let  =  max  |aj^|  and  jh]  =  max  |h,  ) .  Then  we  obtain  the  estimate 
k  k 


(4.6)  |lDF(q)h  -  DF(q*)hllj  <  A^jh|^Ejj  J(t-qj^;q)  -  5(t-q^;q*)  jdt 


n  rtT 


An,lh|  E  r  |y(t-q,  ;q*)  -  y(  t-q,  ;q* )  [dt . 
m  k  k 


Now  from  (4.1)  we  obtain 


T  T 

(4.7)  “  y(t-qj^;q*)  Idt  <  J  |y(t;q)  -  y(t;q*)|dt 

n  „T 

-  ^  ly(t-q.;q)  -  y(t-q^ ;q*) |dt 

“1=^0  J  J 


n  „T 


-  ^  f  |y(t-q.;q)  -  y(t-q. ;q*) |dt 

j=:ro  ••  J 


+  A 


n  T 

^  f  iy(t-q.;q*)  -  y(t-q^;q<')  (dt 
j=l''0  J  J 
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n  T 

<  A  E  f  jy(t;q)  -  y(t;q*)|dt 


n  «T 


+  A 


S  f  |y(t-q.;q*)  -  y{ t-qf ;q* ) | dt . 

j=ro  ^  ^ 


Now  since  y(t;q)  =  S(t;q)xQ  is  differentiable  with  respect  to  q  it  is  not 
difficult  to  show  that  there  are  constants  /?  and  6  such  that 


T 

(4.8)  r  jy(t;q)  -  y(t;q*)|dt  <  f3\q  -  q*] 

■^0 


whenever  jq  -  q* j  <  8.  Combining  (4.7)  and  (4.8)  with  (4.6)  yields 


(4.9)  |jDF(q)h  -  DF(q*)h||^  <  A^|hln^|q  -  q*| 

n  T 

+  A^|h|n  E  f  |y(t-q,  ;q*)  -  y( t-q*;q*) |dt 
n  pT 

+  A  lh|  E  |y(t-q.;q*)  -  y( t-q* ;q*) |dt  . 

”  k=l‘'0  * 


Since  (t?,^)  €  D,  we  have  y  and  y  in  L^(-r,T).  Therefore,  the  integral 

terms  in  (4.9)  approach  zero  as  q  -♦  q*  and  (H8)  holds.  If  ip  is  of  bounded 

variation  on  [-r,0],  then  y  and  y  are  of  bounded  variation  on  [-r,T].  By 
[15,  Theorem  2.1.7(b)]  this  implies  that  the  integral  terms  in  (4.9)  are 
0(  |q  -  q*  1 )  as  q  -»  q*  so  that  (H12)  holds. 

Lemma  4.3.  Suppose  u  €  L^(0,T).  Then  DG(q)  as  defined  by  (4.5)  satisfies 

(HlO).  Moreover,  if  in  addition  u  is  of  bounded  variation  on  [0,T],  then 
DG(q)  satisfies  (H13). 
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Proof :  Using  (4.5)  in  place  of  (4.4)  one  obtains  the  estimate  (4.9)  above 


with  y  replaced  by  z.  Now  if  u  6  L^O.T)  then  z  and  z  are  in  L*(-r,T)  so 
that  (HIO)  holds.  Similarly,  if  u  is  of  bounded  variation  on  [0,TJ,  then  z 

and  z  are  of  bounded  variation  on  [-r.T]  so  that  (H13)  is  satisfied. 

5.  The  Algorithm 

In  this  section  we  define  an  estimation  algorithm  over  a  finite 
dimensional  parameter  space  based  on  quasi  1 ineari zat ion  and  establish  local 
convergence  using  the  results  of  Section  3.  In  particular,  we  assume  that 

the  parameter  space  P  is  P'^with  canonical  basis  e.,  i  =  1,  2,.  .  .,  n. 

This  algorithm  can  also  be  cast  in  a  separable  Hilbert  space  as  in  [17]. 

Given  6  D  and  q  €  P  C  R*'  a  strong  solution  of  (1.1)  is  given  by 

S(t;q)xQ  +  Q(t;q).  Here  we  have  used  the  notation  of  Section  2.  Let  C  be 

I 

a  bounded  linear  mapping  from  X  into  R  and  define 
'lf(t;q)  =  C[S(t;q)xQ  +  Q(t;q)]. 

The  parameter  estimation  algorithm  is  related  to  the  following  optimization 
problem. 

I 

Problem  5.1.  Let  yjeR,j=l.  2,  .  .  .,mbe  data  values  taken  at 
times  t j  6  [0,  T] ,  j  =  1 ,  2,  .  .  . ,  ra,  respectively.  For  q  e  P  define  the 
quadratic  cost  function 

J(q)  =  E  h(  t  -q)  -  y.  1  . 
j=l  •'  •' 

Find  q*  g  P  such  that  J(q*)  <  J(q)  for  all  q  €  P. 
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The  quasi  1 inearization  method  defines  a  recursive  algorithm  whose 
fixed  point  is  a  local  solution  of  Problem  5.1.  A  more  complete 
exposition  is  given  in  [71.  Given  an  initial  guess  q^  €  P  define 

Vi  =  “  -  “.i.a.ii.  -- 

where 

f(q)  =  q  -  [D(q)l  ^b(q) 

D(q)  =  E  M  (t . ;q)M(t . ;q) 
j=l  ^ 

in  « 

b(q)  =  E  M  (t.;q)['y(t  ;q)  -  y  ] 
j=l  J  J  J 

and  the  matrix  M(t;q)  has  its  i^^  column  M^{t;q)  given  by 

M^(t;q)  =  CD^[S(t;q)xQ  +  Q(t;q)]ej,  i  =  1,2,3,. ..,n. 

Lemma  5.1.  Suppose  the  hypotheses  of  Lemmas  3.1  and  3.2  are  satisfied. 
Then  M(tj;q)  is  continuous  in  q  at  q*. 

Proof .  This  is  a  direct  consequence  of  Lemmas  3.1  and  3.2  and  the  above 
definitions . 

Lemma  5.2.  Suppose  the  hypotheses  of  Lemmas  3.3  and  3.4  are  satisfied. 
Then  there  exist  constants  a,  6  >  0  such  that 

jM(t.;q;  -  M(t^;q*)l  <  a|q  -  q*]. 
for  |q  -  q  |  <  <5,  j  =1,2 . m. 
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Proof ■  This  is  a  direct  consequence  of  Lemmas  3.3  and  3.4  and  the  above 


definitions. 

Note  that  although  the  smoothness  results  of  the  previous  sections 
hold  for  an  infinite  dimensional  parameter,  the  implementation  of  the 
solution  of  Problem  5.1  by  this  method  is  limited  to  finitely  many 
parameters.  In  fact  a  simple  rank  argument  is  used  in  [17]  to  show  that 
if  the  number  of  parameters,  n,  exceeds  the  number  of  data  values,  ml. 
then  the  matrix  D(q)  is  singular.  In  [17]  a  pseudo-inverse  is  proposed  as 
a  means  of  solving  the  underdetermined  problem. 

We  can  now  prove  the  following  convergence  results.  These  results  are 
typical  of  quasilinearization  methods  and  the  proofs  given  here  are  in  the 
same  spirit  as  those  in  [7].  We  obtain  superlinear  convergence  when  there 
is  an  exact  fit  to  data  (Theorem  5.1)  and  linear  convergence  in  the 
presence  of  error  (Theorem  5.2). 


Theorem  5,1.  Suppose  the  hypotheses  of  Lemmas  3.1  and  3.2  are  satisfied. 

Moreover,  assume  [D(q*)]  '  exists,  f(q*)  =  q*,  and  J(q*)  =  0.  Then  for 
every  e  >  0  there  exists  S  >  0  such  that 


|f(q)  -  f(q*){  <  £|q  -  q* ] 


for  (q  -  q* (  <  6.  In  particular,  there  is  a  neighborhood  U  of  q*  such  that 
qj^  -♦  q*  as  k  -♦  oe  whenever  q^  e  U. 

Proof.  Note  that  f(q*)  =  q*  implies  that  b(q*)  =  0,  or 

^  rp 

(5.1)  E  'r(t.;q*)['y(t.;q*)  -  y .  ]  =  0. 

j=i  •>  J  J 
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Therefore 


f(q)  -  f(q*)  =  D(q)  ^[3(q)(q  -  q*)  -  b(q)] 


r  ®  — 

=  D(q)"^  E  M^(t.;q)[M(t.:q)(q  -  q*)  -  {'7(t.;q)  -  y.)] 
.j=l  J  •• 


=  D{q)"^  E  M'^(t,;q)[M(t.;q)  -  M(t.;q*)](q  -  q*) 
j=l  J  J  J 


m  _ 

-  D(q)  ^  E  M  (t.;q)[7(t.;q)  -  7(t:;q*)  -  M(t  •q*)(q  -  q*)] 
j-1  J  J  J  J 


-  D(q)~^  E  M^(t .  ;q)[7(t.  ;q*)  -  yj. 
j=l  J  J  J 


Therefore,  using  (5.1)  we  have  that 


(5.2)  f(q)  -  f(q*)  = 


D(q)’^  E  M^(t.;q)[M(t.;q)  -  M(t.;q*)](q  -  q*) 
j=l  ■'  J  J 


-  D(q)''^  E  M^(t.  ;q)[7(t.;q)  -  7(t.;q*)  -  M(t.  ;q)(q  -  q*)] 
j=l  J  J  J  J 


—  1  T  T  - 

-  D(q)  E  [M  (t.;q)  -  M  ( t . ; q* ) ] [7( t . ;q* )  -  y.]. 

j=l  ^  J  J  J 


Note  that  D(q)  *  exists  and  is  bounded  in  a  neighborhood  of  q*  since 
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D(q*)  ^  exists  by  assumption  and  D(q)”^  is  continuous  at  q*  by  Lemma  5.1. 

Let  £  >  0  be  given.  Using  Lemma  5.1  it  is  easy  to  see  that  there 
exist  constants  >0  such  that  the  first  term  in  (5.2)  is  bounded  by 

-  Q*1  for  Iq  -  q*!  <  Furthermore,  since  M(tj;q*)  is  the  Frechet 

derivative  of  7(tj;q)  at  q*,  one  can  show  that  there  exist  constants 

/J  ,  S  >  0  such  that  the  second  term  of  (5.2)  is  bounded  by  f/J  (q  -  q*l  for 

A  *  2 

Iq  -  q* I  <  Combining  these  estimates  with  (5.2)  yields 

(5.3)  |f(q)  -  f(q*)|  < 

f/3jq  -  q*|  +  |D(p)"^|  E  |M^(t.;q)  -  M^(t.;q*)|  h(t.;q*)  -  y J  , 

j=l  •’  J 

for  |q  -  q*|  <  6  =  min  iS  ,6  )  and  0=0.+  0..  Since  J(q*)  =  0,  the  last 

i  2  12 

term  in  (5.3)is  zero.  This  estimate  yields  the  desired  result. 

The  following  theorem  does  not  require  an  exact  fit  to  data,  but  does 
place  some  technical  restrictions  on  the  behaviour  of  M  near  q* .  Note 
that  if  Lemmas  3.3  and  3.4  hold  then  there  exists  2  >  0  such  that  for 
0  <  5  <  5  there  exists  a  constant  K(5)  such  that 

E  iM^'^lt.jq)  -  M'^(t.;q*)|  <  K(^)|q  -  q*| 

j=l  ^  J 

for  |q  -  q*j  <  6.  Let  K*  =  lim  sup  K(^)  and  define 

6  i  0 

(5.4)  A*  =  K*|D(q*)  *(  maxj7(tj;q*)  -  y^ ( . 
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Theorem  5.2.  Suppose  the  hypotheses  of  Lemmas  3.3  and  3.4  are  satisfied. 


Moreover,  assume  lD(q*)]  ^  exists  and  f(q*)  =  q*.  Let  X*  be  defined  by 
(5.4)  and  assume  A*  <  1.  Then  there  exists  6*  >  0  such  that 

|f(q)  -  f(q*)|  <  A*|q  -  q* | 

for  Iq  -  q*l  <  5*.  In  particular,  qj^  as  k  -♦  «  whenever 

IQq  -  q*|  < 

Proof .  This  estimate  is  a  direct  consequence  of  (5.3). 


6.  Numerical  Examples 

In  this  section  we  consider  several  examples  in  which  the  above 
algorithm  was  used  to  solve  parameter  estimation  problems  in  delay- 
differential  equations.  In  these  examples  the  emphasis  is  on  delay 
identification  since  in  the  abstract  setting  this  represents  an  unbounded 
perturbation  of  the  generator  as  noted  in  Section  4. 

With  the  exception  of  Example  6.8,  the  various  unknown  parameters  are 
estimated  using  data  generated  from  closed-form  expressions  for  the 
solution  found  by  the  "method  of  steps".  The  algorithm  is  implemented  by 
an  averaging  scheme  [2]  which  approximates  the  state  equation  and  the 
associated  sensitivity  equations  by  a  system  of  ordinary  differential 
equations.  This  system  is  solved  by  a  fourth-order  Runge-Kutta  routine. 

In  the  one  delay  examples  the  averaging  scheme  is  implemented  with  the 
delay  interval  [-r,0]  divided  into  sixteen  equal  segments,  except  that 
Example  6.8  uses  64  equal  segments.  In  the  two  delay  examples  the 
intervals  [-r2,  -rl]  and  [-rl,0j  are  divided  into  sixteen  equal  segments. 
All  computations  were  done  on  a  VAX  11/750  minicomputer  or  a  SUN 
Microsystem  at  the  Institute  for  Computer  Applications  in  Science  and 
Engineering  (ICASE). 
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Example  6.1.  This  example  illustrates  the  rapid  convergence  of  the  method 

for  a  single  unknown  parameter — the  delay  in  the  following  equation — with 
an  initial  guess  which  is  an  order  of  magnitude  greater  than  the  "true 
value"  of  r  =  1.0.  The  equation  and  the  results  of  the  iteration  are  given 
below. 


x(t)  = 

x(t)  = 

-2x(t)  +  3x(t- 

t  +  1,  t  <  0 

-r),  t  >  1 

iterate 

r 

error 

0 

10.000 

34.056 

1 

1.299 

0.955 

2 

0.946 

0.175 

3 

0.989 

0.115 

4 

0.987 

0.115 

The  convergence  of  the  states  to  ten  data  points  on  the  interval  [0,2]  is 
illustrated  in  Figure  1. 

Elxample  6.2.  The  data  is  the  same  as  for  Example  6.1,  however  in  this  case 

the  algorithm  is  asked  to  estimate  the  coefficients  as  well  as  the  delay. 
The  equation  shows  an  insensitivity  to  the  individual  coefficients  which 
leads  to  the  inaccuracy  in  the  converged  estimates.  In  fact,  because  of 
errors  introduced  by  the  averaging  scheme  for  computing  the  state,  the 
estimated  values  fit  the  data  better  than  the  "true  values"  used  to  compute 
the  data  by  the  method  of  steps.  The  "true  values"  are  a  =  -2,  b  =  3,  and 
r  =  1.  The  equation  and  the  results  of  the  iteration  are  given  below: 


f  x(t)  =  ax(t)  +  bx{t-r),  t  >  0 
I  x(t)  =  t  +  1,  t  <  0 
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iterate 

a 

b 

r 

error 

0 

-4.000 

7.000 

2.000 

3.379 

1 

-0.815 

3.537 

1.184 

2.968 

2 

-1.596 

3.342 

1.122 

0.775 

3 

-2.403 

3.713 

1.002 

0.188 

4 

-2.250 

3.361 

1.015 

0.094 

5 

-2.352 

3.483 

1.006 

0.093 

The  convergence  of  the  states  is  illustrated  in  Figure  2. 

Example  6.3.  This  case  illustrates  the  effect  of  a  forcing  function  on 
state  equation.  The  nonhomogeneous  delay-differential  equation 


I  x(t)  =  ax(t)  +  bx(t-r)  +  u(t),  t  >  0 
^  x(t)  =  t  +  1,  t  <  0 


where 


u(  t) 


'  0  ,  t  <  0.1 
1  ,  t  >  0.1 


is  solved  in  closed  form  by  the  method  of  steps  with  parameter  values 
a  -  -2,  b  =  3,  r  =  1  as  in  Example  6.2.  The  results  of  the  parameter 
estimation  algorithm  are  given  below: 


iterate 

a 

b 

r 

error 

0 

-4.000 

7.000 

2.000 

4.0527 

1 

1.022 

3.165 

1.140 

39.2657 

2 

-2.637 

23.652 

1.168 

24.9577 

3 

-5.979 

28.631 

1.141 

11.6964 

4 

-8.034 

23.250 

1.118 

3 . 5425 

5 

-5.167 

5.417 

1.028 

2.0471 

6 

-1.239 

4.195 

1.008 

4.8981 

7 

-2.861 

6.222 

1.005 

1 . 8930 

8 

-2.485 

3.795 

0.998 

0.0819 

9 

-2.115 

3.201 

1.013 

0.0724 

10 

-2.247 

3.380 

0.998 

0.0691 

The  results  are  similar  to  those  of  Example  6.3,  except  that  the  solution 
has  become  somewhat  more  sensitive  to  the  coefficients. 

Example  6.4.  This  example  indicates  the  ability  of  the  algorithm  to 

estimate  two  unknown  delays,  llje  algorithm  converges  rapidly  from  a 
relatively  poor  initial  guess.  The  "true  values”  are  r^  =  1.0  and 

r^  =  2.0.  The  equation  and  the  results  of  the  parameter  estimation 

algorithm  are  given  below  and  the  convergence  of  the  states  to  ten  data 
points  on  the  interval  [0,3]  is  illustrated  in  Figure  3. 


x(t)  =  -x(t)  +  x(t-rj)  -  x(t-r2),  t  >  0 
x(t)=t+l,  t<0 


iterate 

•■l 

^2 

error 

0 

0.600 

4.000 

7.500 

1 

1.569 

3.216 

2.295 

2 

1.146 

2.100 

0.100 

3 

0.977 

1.998 

0.034 

4 

0.978 

2.003 

0.032 
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Example  6.5.  The  equation  and  data  for  this  example  are  the  same  as  in 

Example  6.4.  In  this  case  the  initial  guess  reverses  the  order  of  the 
"true"  delay  values.  The  results  of  this  iteration  are  given  below  and 
covergence  of  the  states  on  the  interval  [0,3]  is  illustrated  in  Figure  4. 


i terate 

*■2 

error 

0 

2.000 

1.000 

2.460 

1 

0.483 

1.151 

1.379 

2 

1.561 

2.014 

0.788 

3 

1.100 

2.072 

0.077 

4 

0.980 

2.002 

0.033 

Example  6.6.  In  this  case  the  algorithm  is  asked  to  estimate  parameters  in 

a  delay  model  of  a  system  with  no  delay.  Ten  data  points  on  the  interval 
[0,2]  are  computed  from  the  exponential  solution  of 

/  x(t)  =  -•2x(t) 

\  x(0)  =  1 


and  the  algorithm  is  asked  to  estimate  unknown  parameters  in  the  system 

f  x(t)  =  ax(t)  +  bx(t-r),  t  >  0 
I  x(t)  =  t  +  1,  t  <  0 

The  first  four  iterations  are  given  below; 


i terate 

a 

b 

r 

error 

0 

-3.000 

3.000 

2.000 

1.2577 

1 

-3.060 

-0.637 

1.947 

0.2551 

2 

-1.687 

0.235 

1.981 

0.1144 

3 

-1.967 

0.025 

1.985 

0.0110 

4 

-2.000 

0.000 

1.986 

0.0001 

On  the  fifth  iteration  the  algorithm  aborted  when  it  was  asked  to  invert  a 
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nearly  singular  matrix.  'Hiis  reflects  the  fact  that  at  the  true  parameter 
values  the  state  is  completely  insensitive  to  the  delay. 

Example  6.7.  This  case  is  the  same  as  the  previous  example  except  that  the 

data  is  taken  from  the  closed  form  solution  of  the  nonhomogeneous  undelayed 
equation 

/  x(t)  =  -2x(t)  +  u(t) 

I  x(0)  =  1 

where  u  is  the  same  step  function  as  in  Example  6.3.  The  results  are 
similar  to  those  of  the  previous  example. 


iterate 

a 

b 

r 

error 

0 

-3.000 

3.000 

2.000 

1.3135 

1 

-2.848 

0.099 

1.804 

0.5121 

2 

-1.841 

0.138 

2.401 

0.0811 

3 

-1.971 

0.003 

2.508 

0.0197 

Example  6.8.  In  this  example  we  consider  the  second-order  equation 

^-|■{t)  +  w^x(t)  +  a^  ^U-r)  +  ajx(t-r)  =  u(t),  t  >  0, 
dt 

,  x(t)  =1,  t  <  0, 

where  u(t)  is  the  step  function  of  Example  6.3.  This  equation  models  a 
harmonic  oscillator  with  retarded  damping  and  restoring  forces.  In  [13]  a 
quasi  linearization  algorithm  is  used  to  estimate  coefficients  in  this 
equation.  The  methods  of  this  paper  allow  the  delay  r  to  be  added  to  the 
set  of  unknown  parameters.  For  this  example  the  averaging  method  was  used 
to  compute  "data"  values  for  the  parameter  estimation  algorithm  with  "true" 
values  of  w  =  6,  Sq  =  2.5,  a^  =  9,  and  r  =  1.  The  results  of  the  iterative 

algorithm  are  given  below  and  the  convergence  of  the  states  (displacement 
and  velocity)  on  the  interval  [0,  2]  is  illustrated  in  Figures  5  and  6. 
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iterate 

u 

®0 

“l 

r 

error 

0 

4.100 

4.600 

6.300 

1.500 

15.212 

1 

5.073 

6.025 

-8.338 

0.918 

15.181 

2 

6.705 

4.710 

-0 . 682 

1.524 

12.389 

3 

6.188 

-14.677 

-4.838 

1.102 

31.950 

4 

5.902 

12.347 

8.396 

1.068 

25.234 

5 

5.964 

2.994 

8.980 

1.061 

2.186 

6 

5.995 

2.416 

9.016 

1.004 

0.344 

7 

6.000 

2.503 

8.999 

1.000 

0.007 
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Abstract.  We  consider  identification  of  parameters  in  a 
Volterra  integro-differential  system  with  a  weakly  singular 
kernel.  Such  kernels  arise  in  fractional  derivative  damping 
models  of  viscoelastic  materials.  The  Volterra  equation  is 
cast  in  a  semigroup  setting  to  establish  results  on  the 
differentiability  of  the  solution  with  respect  to  a  parameter. 
These  results  are  needed  for  convergence  of  the  identification 
algorithm.  Numerical  results  are  presented. 


1.  Introduction.  In  this  paper  we  consider  the  identification 
of  parameters  in  a  Volterra  integro-differential  equation  with 
a  singular  kernel.  The  equation  of  interest  has  the  form 


(1.1) 


r  ^ 

'  w(t)  =  Mw(t)  +  I  K(t-s,p)w(s)ds  +  F(t),  t  £  0, 

<  ••CO 

^  W(0)  =  7J,  W(s)  =  <^(S),  S  <  0, 


^This  research  was  supported  by  the  Air  Force  Office  of 
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where  M  is  an  n  x  n  constant  matrix,  17  €  r”,  0  €  L^(-<»,0;IR*^) 
and  K(*,p)  is  an  n  x  n  singular  kernel  depending  on  a 
parameter  p  contained  in  an  admissible  parameter  set.  We  are 
particularly  interested  in  a  kernel  function  of  the  form 


g(s.p) 


r(l-a)s“ 


s  >  0, 


.  3 

where  r(*)  denotes  the  gamma  function,  p  =  (o£,0,r)  e  IR  with 
0  s  a  <  1  and  /S  >  0.  Such  kernels  arise  in  the  study  of 
fractional  derivative  models  of  viscoelastic  structures.  For  a 
more  complete  discussion  of  the  origins  of  this  kernel  and  the 
viscoelastic  models  we  refer  the  reader  to  [12],  [18],  [13], 
[15],  and  in  particular  to  [17]  and  the  extensive  bibliography 
therein. 

Banks,  et.al.  [2]  have  identified  parameters 
corresponding  to  ^  and  r  in  a  similar  (but  different)  model, 
but  assumed  that  a  was  known.  Torvik  and  Bagley  [1],  [18]  have 
estimated  the  parameter  a,  but  in  the  Laplace  transform 
domain.  In  this  paper  we  restrict  ourselves  to  identifying  a 
only,  though  the  theory  may  be  modified  to  include  B  and  y  as 
well. 

In  order  to  relate  equation  (1.1)  to  a  (idealized) 
physical  model,  consider  the  longitudal  motions  of  a  uniform 
bar  fixed  at  both  ends  with  Boltzmann  type  damping.  The 
governing  equation  is  ([9],  [14]) 


(1.2) 


f  pu^^(x,t) 


—  |EUj^(x,t)  +  J^g(t-s)u^(x,s)ds 


} 


+  f(X,t),  0  <  X  <  1,  t  >  0, 


with  boundary  conditions 
and  initial  conditions 


u(0,t)  =  0, 
u(x,0)  =  d(x) , 


u(l,t)  =  0, 
u^(x,0)  =  v(x) . 


Here,  u(x,t)  represents  the  axial  displacement  of  position  x 
at  time  t,  p  is  the  density  of  the  material,  E  a  stiffness 
parameter,  f(x,t)  a  forcing  function,  and 


g(s) 


r(l-a)s“ 


represents  a  fractional  derivative  damping  term  modified  to 
have  exponential  decay  [12]. 

A  common  approach  to  the  parameter  identification  problem 
[4]  is  to  apply  a  Galerkin-type  approximation  scheme  to  the 
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beam  equation  and  then  incorporate  some  type  of  identification 
algorithm  to  the  approximating  system  of  integro-dif ferential 
equations.  If  one  applies  a  Galerkin  scheme  to  equation  (1.2) 
(e.g.  using  linear  splines) ,  one  obtains  a  system  of  equations 
of  the  form 


(1.3) 


-  A^-^(t)  =  Bv(t)  +  f  g(t-s)v(s)ds  +  F(t)  . 

dt‘‘  ■'o 


In  this  equation  A,  B,  and  C  represent  constant  matrices  and 
v(t)  and  F(t)  are  vectors  of  appropriate  dimension. 

In  order  to  retain  the  salient  features  but  simplify  the 
analysis  in  the  following  sections,  we  shall  consider  the 
following  scalar  version  of  (1.3)  : 


(1.4) 


f  d^x, 

dt‘ 


(t)  =  ax(t)  +  Sr  f  g(t-s)x(s)ds  +  f(t), 

J 


[  x(0)  =  Xq,  x(0)  = 


^1  • 


Integrating  (1.4)  we  obtain 
■ 

x(t)  ==  a  x(s)ds  +  g(t-s)x(s)ds  +  f(t), 


(1.5) 


where 


x(0)  =  Xq, 


f(t)  =  Xj^  +  J  f(s)ds 


Define  z(t)  =  x(s)ds  ,  then  z(t)  =  x(t)  and  we  obtain 

J  0 

the  system  of  integro-dif ferential  equations 


(1.6) 


x(t) 

2(t) 


=  az(t)  +  g(t-s)x{s)ds  +  f(t), 
0 

=  x(t)  , 


with  x(0)  =  Xq,  z(0)  =  0. 

A  standard  assumption  in  viscoelasticity  [13]  is  that  the 
material  is  in  an  unstrained  state  for  time  t  <  0.  This  would 
correspond  to  u(x,s)  =  0  for  s  <  0  in  equation  (1.2).  It 
follows  then  that  x(s)  =  0  and  2(s)  =  0  for  s  <  0  in  (1.6).  If 
we  define  w(t)  =  col (x(t) , z (t) ) ,  then  w(t)  satisfies 


5 


(1.7) 


w(t)  =  Mw{t)  +  K(t-s)w(s)ds  +  F(t), 
''O 

w(0)  = 


where  M  =  (  J  q  )'  “  (  ^0^^  o  )  '  “  (  0 

Since  w(s)  =  col (0,0)  for  s  <  0,  we  may  rewrite  (1.7)  as 


f(t) 


I 


w(t)  =  Mw(t)  +  K(t-s)w(s)ds  +  F(t) 


w(0) 


s  <  0, 


which  is  in  the  form  of  equation  (1.1). 

The  remainder  of  the  paper  is  outlined  as  follows.  In 
Section  2  we  review  previous  results  that  place  equation  (1.1) 
in  a  semigroup  setting  in  order  to  establish  existence  of 
solutions.  Differentiability  results  needed  for  the  parameter 
estimation  algorithm  are  then  proved.  In  Section  3  the 
quasilinearization  algorithm  used  for  the  identification 
procedure  is  discussed  along  with  convergence  results. 
Numerical  examples  are  given  in  Section  4 . 


2.  The  abstract  setting.  In  this  section  we  develop  an 
abstract  framework  for  the  Volterra  integral  equation 
discussed  in  the  previous  section.  Namely,  we  will  consider 
equation  (1.1)  in  the  form 


(2.1) 


where  t)  = 


'  w(t) 

=  Mw(t)  + 

r  K(-s,a)w(t+s)ds  +  F(t),  t  >  0, 

J  -00 

w(0) 

=  V,  W(s) 

=  (p(S)  ,  s  <  0, 

X 

.0  . 

€  IR^,  M  = 

( ?  ? )'  '  (  T’]' 

,  s  >  0. 
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We  assume  ^  is  a  positive  constant  and  0  s  a  <  i.  By  a 
solution  of  (2.1)  we  mean  a  function  w:  (-»,«)  -»  such  that 

w  is  absolutely  continuous  (A.C.)  on  [O,®)  and  satisfies  the 
integral  equation  a.e.  on  /  w(0)  =  tj,  and  w(s)  =  ip(s) 

a.e.  on  (-<*>,0]. 

Our  semigroup  formulation  follows  the  construction  in  [5] 
as  further  developed  in  [10]  and  [11]-  Define  the  product 
space  X  =  K  X  L^(-oo,o)  with  norm 

•  Consider  the  homogeneous 

A  •"  \  '“®  <  ^  / 

equation 


f  y(t)  =  My(t)  +  K(-s,a)y (t+s)ds,  t  >  0, 

(2.3)  j 

{  y(o)  =  V,  y(s)  =  ip(s) ,  s  <  0. 

Then  for  each  pair  (v,<p)  e  X,  (2.3)  has  a  unique  solution,  and 
moreover  the  mapping  S(t,a)(T),9)  =  (y(t),y^(-))  defines  a 

CQ-semigroup  on  X.  Here  we  have  used  the  notation 

=  y(t+s)  ,  t  £  0,  s  <  0. 

Fix  c  e  (0,  1)  and  define  the  parameter  set  P  =  [0,  1-e]. 
Then  it  is  readily  seen  from  (2.2)  that  there  is  a  constant  C, 
independent  of  a,  such  that 

09 

(2.4)  r  |K(s,a)lds  £  C  for  all  a  e  P. 

0 

Under  this  condition  it  is  shown  in  [10]  and  [11]  that  the 
semigroup  S(t,a)  is  generated  by  a  closed  and  densely-defined 
operator  i4(a)  defined  by 

Dom(j4(a))  2  D  =  [{T\,(p)  €  X;  (p  is  A.C.  on  compact  subsets 

of  (-co,0],  (p  e  L*(-a),0),  ^(0)  =  TJ) 

and 

fO 

«*(«)('r?,^)  =  (Mtj  +  K{-s,a)(p(s)ds,  <p  )  for  e  D. 

''  -CO 

Our  task  is  to  show  that  the  solution  w(t,a)  s  w(t)  of 
(2.1)  is  differentiable  with  respect  to  a  and  that  this 
derivative  is  sufficiently  smooth  to  establish  the  local 
convergence  of  the  algorithm  defined  in  Section  3.  This 
involves  verifying  the  conditions  in  the  semigroup  setting 
established  in  [6]  and  [8].  We  therefore  assume  in  what 
follows  that  the  reader  has  these  papers  in  hand. 

Since  we  are  interested  in  dependence  on  a,  we  write 


jtf(a)  =  +  S(a)  where  d  is  independent  of  a  and 

(2.5)  ^{a){-n,(p)  =  (  I  K(-s,a)^(s)ds,  0  ),  €  D. 

—00 

Note  that  the  range  of  S(a)  is  the  finite-dimensional  space 
Y  =  R  X  {  0 )  . 

Fix  Yq  €  X,  €  P,  and  T  >  0.  Then  the 
differentiability  with  respect  to  a  at  of  the  solution 
y(t,a)  of  (2.3)  is  a  consequence  of  the  following  theorem. 

Theorem  2.1.  For  every  t  e  [0,T3,  S(t,a)yQ  as  defined  above 
is  Frechet  differentiable  with  respect  to  a  at  and  its 
derivative  is  given  by 

D^S(t,ao)yo  "  J^S(t-s,aQ) [D^5(aQ)yQ] (s)ds,  0  s  t  s  T, 

where  ^(a)  is  for  each  a  e  P  a  mapping  from  X  into  L^(0,T;X) 
defined  by 

(2.6)  [^(a)yo](t)  =  (  I  K(-s,a)y(t+s)ds,  0  ) 

—00 

for  Yq  e  X,  0  s  t  s  T.  Recall  that  y  is  the  solution  of  (2.3) 
with  (Tj,?))  =  Yg. 

Proof;  This  result  is  proved  in  [7]  for  a  general  Volterra 
kernel  K(s,a)  satisfying  condition  (2.4)  under  the  following 
hypothesis; 

the  mapping  a  K(*,a)  from  P  into  L^{0,oo)  is  Frechet 
'  ‘  '  differentiable  with  respect  a  at  a^. 

Recall  that  K(s,a)  =  g(s,a) f  ^  ^  j  where 


g(s,a)  =  — ! - -  ,  s  >  0,  a  e  P. 

r(l-a)s“ 

Let  '  denote  differentiation  with  respect  to  a.  Then 
computation  shows  that  g'  and  g"  are  of  the  form 

g'(s,a)  =  g^(a)(ln  s)e'^®s‘“  +  g2{a)e"^®s‘“ 


48 


and 


g"(s,a)  =  g2(a)(ln  +  g^(a)(ln  s)e  -"s 

+  g^CoOe-^^s-o 


-gs-a 


Where  g^^, .  .  .,gg  are  continuous  functions  of  a  on  P  which  can 

be  exp] icitly  calculated  in  terms  of  the  gamma  function  and 
its  derivatives.  Important  properties  of  g'  and  g"  for  our 
purposes  are  that  there  are  functions  ijt.,  if/  c  L^{0,oo)  such 

that 

(2.8)  Ig' (s,a)  I  s  for  s  >  0,  a  e  P, 

and 

(2.9)  lg"(s,a)|  s  i/'2(s)  for  s  >  0,  a  €  P. 


For  example,  one  can  take 

'  (Mjln  sle"^®  4-  M_e"^®)/s^”^,  for  0  <  s  <  1 

=  ■  -fis  -fls 

[  Mj^lln  s|e  ^  +  M2e  for  sal, 

where  and  M2  are  upper  bounds  on  P  of  |g^(a)|  and  |g2(a)|, 
respectively.  There  is  an  analogous  expression  for 
Therefore,  by  Taylor's  theorem  with  remainder,  we  obtain 


|K(s,a  +  h)  -  K(s,a)  -  K'(s,a)hl 

=  lg(s,a  +  h)  -  g(s,a)  -  g- (s,a)h| 

=  lg"(s,Ci(s))h2/2| 
s  111^(5)  IhiVz 

for  s>0,  a,  a  +  heP,  and  Cj^(s)  between  a  and  a  +  h. 

Integrating  this  inequality  over  (0,<»)  and  using  e  L^(0,oo) 

yields  (2.7)  and  completes  the  proof  of  Theorem  2.1. 

A  sufficient  smoothness  property  for  the  local 
convergence  of  the  parameter  estimation  algorithm  is 
established  in  the  following  theorem. 


Theorem  2.2.  For  every  t  €  [0,T3,  a*  €  P  and  y^  e  x, 
D^S(t,a)yQ  is  strongly  locally  Lipschitz  continuous  with 
respect  to  a  at  a*. 
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Proof:  The  proof  relies  on  Lemma  3.3  of  [8].  We  must  show 
that  hypotheses  (Hll)  and  (H12)  of  that  paper  hold  in  this 
application.  By  definition  (2.6),  (Hll)  requires  that  there 
exist  constants  5^  >  0  such  that 


T  0 

(2.10)  I  ||  (K(-s,o£+h)  -  K(-s,a)  )y  (t+s)  dsidt  s  K^IhllTjl 

0  ”00 


for  Ihl  s  5^,  where  y  is  the  solution  of  (2.3)  with  a  =  a*  and 

(p  s  0.  Note  that  since  (y(t)  ,y.)  =  S  (t ,«. )  (tj,  0)  we  have 

|y(t)  I  s  M^e  iTjl  for  t  £  0,  and  y(t)  =  ^>(t)  =  0  for  t  <  0. 

It  is  shown  in  (7]  that  the  constants  and  o)  may  be  taken 

independently  of  a  e  P.  Therefore,  by  Fubini's  theorem  and 
the  mean  value  theorem  we  obtain 


J  r  |(K(-s,a+h)  -  K(-s,a) )y(t+s) Ids  dt 
0  ‘"00 

fO  T+s 

=  J  |K(-s,a+h)  -  K(-s,a) |  j  ly(t) |dt  ds 

•oo  S 

fO  T 

s  lK(-s,a+h)  -  K(-s,a)|  |y(t) |dt  ds 

•'  0 


s  M^Te“'^|7?|  J  |K(-s,a+h)  -  K(-s,a)  Ids 


s  M^Te‘‘^'’^|T7!lhl  f  ^^(-s)ds 

for  a,  a  +  h  €  P.  Here  we  have  used  (2.8)  in  the  last 
inequality.  Since  e  L^(0,oo),  this  establishes  (2.10). 

Hypothesis  (H12)  of  [8]  requires  the  Lipschitz  continuity 
of  the  derivative  with  respect  to  a  of  the  mappinq  ?(a) 
defined  by  (2.6).  For  brevity,  we  denote  the  value  of  this 
derivative  at  a  e  P  by  D?(a).  The  existence  of  D?(a)  was 
shown  in  Theorem  2.1  and  from  the  proof  of  that  theorem  we 
have  the  formula 

[D?(a)h](t)  =  (  J  [K' (-s,a)h]y(t+s)ds,  0  ) 

”09 

for  0  s  t  s  T,  h  e  R,  a  €  P,  where  y  is  the  solution  of  (2.3) 

I 
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with  a  =  and  (v,9>)  -  Yq'  Recall  that  '  denotes 

differentiation  with  respect  to  a.  The  local  Lipschitz 
continuity  of  D^(a)  at  a  point  a  =  a*  e  P  now  follows  from 
estimates  similar  to  those  used  to  establish  (Hll)  but  with 

in  place  of  This  completes  the  proof  of  Theorem  2.2. 

We  now  turn  our  attention  to  solutions  of  the 
nonhomogeneous  equation  (2.1).  It  is  well-known  that  a  mild 
solution  to  this  equation  is  given  by  the  variation  of 
constants  formula 

(w(t),w^)  =  S(t, a)  (■n,(p)  +  Q(t,a) 

where 

(2.11)  Q(t,a)  =  S(t-s,a) (F(s) ,  0)ds. 

''  0 

It  remains,  therefore,  to  consider  the  existence  and 
smoothness  of  the  derivative  of  Q(t,a)  with  respect  to  a.  We 
again  appeal  to  [8]  where  these  properties  of  Q(t,a)  are 
demonstrated  by  considering  similar  properties  of  the  mapping 
^(a):  L^(0,T;X)  -*  l'(0,T;X)  defined  by 

[&(a)v](t)  =  J^«(a)S(t-s,aQ)v(s)ds 

for  V  €  L^(0,T;X),  a  €  P,  fixed.  Note  that  if 
v(t)  =  (F(t),  0)  for  some  Fe  L*(0,T),  then 

(w(t),w^)  =  J^S(t-s,aQ)v(s)ds 

where  w  is  a  mild  solution  of  (2.1)  with  (Vf<p)  =  (0,0).  Since 
$(a)  is  a  difference  of  closed  operators, 

[5(a)v](t)  =  S(a) (w(t) ,w^) . 

Therefore,  using  definition  (2.5)  we  obtain  in  this  setting 
that 

rO 

(2.12)  [^(a)v](t)  =  (  I  K(-s,a)w(t+s)ds,  0  ) 

—00 

where  v{t)  =  (F(t),  0)  and  w  is  the  solution  of  (2.1)  with 
a  =  Uq  and  (v,<p)  =  (0,0).  Here  we  assume  F  is  sufficiently 


smooth  for  the  solution  w  to  exist  in  the  strong  sense  defined 
earlier. 

Comparing  definitions  (2.6)  and  (2.12),  we  see  that 
properties  of  5(a)  with  respect  to  a  can  be  proven  in  the  same 
way  as  the  corresponding  properties  of  ?(a)  using  the  solution 
of  (2.1)  in  place  of  the  solution  of  (2.3)  at  a  fixed  a^  €  P. 

For  this  reason  the  following  theorems,  which  are  consequences 
of  Lemmas  3.2  and  3.4  of  [8],  are  stated  without  proof. 

Theorem  2.3.  For  F  e  L^(0,as),  let  Q(t,a)  be  defined  by 
(2.11).  Then  for  every  t  e  [0,T]  and  a^  e  P,  Q(t,a)  is 

Frechet  differentiable  with  respect  to  a  at  and  this 

derivative  is  given  by  the  formula 

D^Q(t,ao)  =  J^S(t-s,aQ) [D^5(aQ)v3 (s)ds 

where  v(t)  =  (F(t) ,  0)  and  5(a)  is  defined  by  (2.12). 

Theorem  2.4.  Suppose  the  hypotheses  of  Theorem  2.3  hold. 

Then  the  mapping  D^Q(t,a)  is  locally  Lipschitz  continuous  with 

respect  to  a  at  every  e  P. 

3,  The  algorithm.  In  this  section  we  define  a  parameter 
estimation  algorithm  based  on  quasilinearization  and  state 
some  local  convergence  results.  For  later  adaptation  we 
develop  the  algorithm  for  the  case  a  e  P  c  r"  with  canonical 
basis  e^^,  i  =  1,2,.  .  .,n.  In  Section  4  the  algorithm  is 

applied  in  the  case  n  =  1.  The  definitions  and  theorems 
stated  here  may  also  be  found  in  [8],  but  are  included  here 
for  completeness. 

Using  the  notation  of  the  previous  section,  let 

Yn  “  (V,<p)  €  X  and  a  e  P.  Let  C  be  a  bounded  linear  mapping 

^  t 

from  X  into  a  f inite-dimensional  space  R  ,  and  define 

w(t,a)  =  C[S(t,a)yQ  +  Q(t,a)]. 


The  parameter  estimation  algorithm  is  related  to  the  following 
optimization  problem. 

p 

Problem  3.1.  Let  w^eR,  j=l,2,,  .  .,mbe  data  values 
taken  at  times  tj  e  [0,T],  j  =  1,2,.  .  .  ,m,  respectively.  For 
a  e  P  define  the  quadratic  cost  function 


m 

J(a)  =  Y  lw(t^,a)  >  Wjl^. 
j=l 


Find  a*  €  P  such  that  J(a*)  s  J(o£)  for  all  a  e  P. 


The  quasilinearization  algorithm  method  defines  a 
recursive  algorithm  whose  fixed  point  is  a  local  solution  of 
Problem  3.1.  A  more  complete  exposition  is  given  in  [3]. 
Given  an  initial  guess  €  P  define 

“k+1  ^  ,  k  =  0,1,2, .  .  . 

where 

Ha)  =  a  -  [D(a)  ]"'b(a) 
m 

D(a)  =  Y  M^(tj,a)M(tj,a) 

j=l 


m 

b(a)  =  ^  M^(tj,a)  [w(tj,a)  -  w^] 
j=i 

and  the  matrix  M(t,a)  has  its  i*^*'  column  M^(t,a)  given  by 

M^(t,a)  =  CD^[S(t,a)yQ  +  Q(t,a)]€^,  i  =  1,2,.  .  .n. 

The  following  theorems  are  typical  of  quasilinearization 
methods.  Their  proofs  may  be  found  in  [8].  We  obtain 
super linear  convergence  when  there  is  an  exact  fit  to  data 
(Theorem  3.1)  and  linear  convergence  in  the  presence  of  error 
(Theorem  3.2) . 

Theorem  3.1.  Suppose  the  conditions  of^the  previous  section 
are  satisfied.  Moreover,  assume  [D(a)]  exists,  ({a*)  =  <x* , 
and  J(a*)  =  0.  Then  for  every  e  >  0  there  exists  6  >  0  such 
that 


\Ha)  -  t{a*)\  s  cla  -  a*| 

for  la  -  a*l  s  5.  in  particular,  there  is  a  neighborhood  ll  of 
a*  such  that  aj^  ^  a*  as  k  ->  <»  whenever  e  I/. 

The  following  theorem  does  not  require  an  exact  fit  to 
data  but  does  place  some  technical  restrictions  on  the 
behaviour  of  the  matrix  M(t,a)  near  a*.  Note  that  under  the 
conditions  of  Theorem  2.2  there  exists  number  6  >  0  such  that 
for  0  <  6  <  6  there  exists  a  constant  K(6)  such  that 
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m 

^  |M^(tj,a)  ~  M'^(tj,a*)l  s  K(5)|a  -  a*! 

j=l 

for  la  -  a*l  s  6.  Let  K*  =  lim  sup  K(5)  and  define 

3  0 

X*  =  K*lD(a*)”*|  n»ax|w(t .  ,a*)  -  w .  1 . 

j  ^  ^ 

Theorem  3.2.  Suppose  the  conditions  of  the  previous  section 
are  satisfied.  Moreover,  assume  [D(a*)]~^  exists  and 
f(a*)  =  a*.  Let  X*  be  defined  as  above  and  assume  x*  <  1. 
Then  there  exists  3*  >  0  such  that 

|^(a)  -  f(a*)  I  s  X*la  -  a*| 

for  la  -  a* I  s  3*.  In  particular,  aj^  a*  as  k  «  whenever 
la^  -  a*|  s  5*. 


4.  Numerical  results.  In  this  section  we  present  several 
examples  that  illustrate  the  ideas  discussed  in  the  previous 
sections.  Recall  the  identification  problem  :  given 
observations  at  times  t^,  j  *  l,2,...,m,  determine 

a  €  tO/1)  that  minimizes  tlj^  cost  functional 


J(a) 


where  x(t)  satisfies 


I  (X(tj)  -  w.)2 

j=l 


(4.1) 


X(t) 

X(0) 


+ 


y 

r{i-a) 


t  -0(t-s) 

- —  x(s)ds  +  f(t), 

0  (t~s) 


X 


0  • 


The  quasilinearization  algorithm  requires  that  we  solve  (4.1) 
along  with  its  sensitivity  equation  which  has  the  form 
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°  +  rrW  I 


(4.2) 


t  g-0(t-s) 


+  yr' (i~a)  r 

r(l-a)^  "^0 


t  g>^(t-s) 


(t-s) 


0  (t-s) 


x(s}ds 


a  a 


x^(s)ds 


■  r(l-a)  J 


^  ln(t-s)e-^<^-^^ 


(t-s) 


a 


x(s)ds. 


[  x^(0)  =  0, 

where  x  (t)  =  l^(t)  .  The  zero  initial  condition  reflects  the 
a  dOL 

fact  that  the  value  x(0)  =  is  independent  of  the  parameter 


ct. 

The  implementaton  of  the  identification  scheme  begins 
with  an  initial  guess  for  a.  Equations  (4.1)  and  (4.2)  are 
integrated  using  this  initial  value,  then  x(t)  and  x^{t)  are 

used  to  give  an  updated  estimate  of  the  parameter.  For  this 
particular  problem  the  quasilinearization  algorithm  updates 
the  current  estimate  aj^  according  to 


m 


“)C41  '  “iC  ' 


m 


where  x(t)  and  x^(t)  are  the  solutions  of  (4.1)-(4.2)  computed 

for  a  =  In  order  to  numerically  integrate  the  state  and 

sensitivity  equations  we  first  convert  (4.1) -(4. 2)  to  integral 
equations  via  the  substitution  z(t)  =  x(t)  and  Zj^(t)  =  x^(t). 

Then  one  has  the  4x4  system  of  integral  equations  consisting 
of  (4.1) -(4!,  2)  with  the  above  substitutions  coupled  with 
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where  tp  e  [0,t)  is  determined  by  the  approximation  scheme. 

The  solution  of  the  system  of  the  integral  equations  is  then 
approximated  by  applying  a  product  integration  method  based  on 
Simpson's  rule  to  the  singular  integral  terms,  and  Simpson's 
rule  to  all  other  integrals.  For  a  description  of  product 
integration  methods  we  refer  the  interested  reader  to  [16]. 

In  each  of  the  following  examples  we  numerically  solve 

(4.1) -(4.3)  in  time  on  the  interval  [0,1].  Define  t^  =  j/N, 

j  =  0,1, ...,N.  The  numerical  integration  scheme  then  computes 
values  for  x(tj).  Examples  (4.1)  and  (4.2)  presented  here  are 

computed  using  a  value  of  N  =  50,  and  Example  (4.3)  is 
computed  using  N  =  200.  In  each  case  5  data  points  located  at 
t=.2,  .4,  .6,  .8,  and  1.0  are  used  in  the  identification 
procedure.  The  true  values  of  x(t)  in  all  of  the  Figures 

(4.1)  “(4. 5)  are  denoted  by  x. 


Example  4.1.  In  this  example  we  set  the  values  of  a,  $,  and  r 
to  1.,  1.,  and  5. respectively.  The  parameter  value  to  be 
identified  is  a  =  ^  ,  and  the  nonhomogeneous  term  f(t)  is 


f(t)  =  e“^(  1024t^  +  2176t  +  1792  )  -  1825 

(  32768  .4  8192  .3  512  ^2 

r(.5)  [  315  ^  ”  35  3 


128 

3 


t  +  2 


.  .  — t 

The  true  solution  is  x(t)  =  e  T^(2t  -  1)  where  T^(s)  is  the 

Chebyshev  polynomial  of  degree  4  on  -1  s  s  s  i.  Tables  (4.1) 
and  (4.2)  contain  the  results  for  two  computer  runs,  one  for 
an  initial  a  of  =  .9,  and  the  other  for  =  .2  .  The 

sequence  of  values,  their  corresponding  costs  J(aj^),  and 

the  values  of  the  state  x(t)  at  time  t  =  1  are  included  to 
illustrate  the  convergence.  The  true  value  of  x(l)  is 
.3678794.  Note  that  in  each  case  once  an  estimate  of  a  is 
greater  than  .5,  then  the  sequence  of  iterates  converges 
monotonically  down  to  the  true  value.  This  is  a  characteristic 
of  all  of  our  computer  simulations,  and  seems  to  indicate  that 
it  is  better  to  choose  an  initial  value  of  a  that  is  high 
instead  of  low.  In  fact,  for  all  simulations  another 
characteristic  is  that  if  the  initial  choice  is  excessively 
low,  then  the  next  estimate  of  a  is  greater  than  1,  and  the 
integral  becomes  undefined.  For  this  particular  example, 
aQ  =  . 1  resulted  in  a  value  greater  than  1  for  the  next 

iterate  and  the  integration  scheme  broke  down.  However,  though 
not  shown  here,  some  examples  ran  succesfully  even  with  a 
negative  initial  value  for  a.  Figures  (4.1) -(4. 2)  show  the 
convergence  of  the  state  x(t)  for  the  initial  values  of 
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Oq  =  .9,  and  .2,  respectively. 


5,  0=1,  true  a 

=  .5,  initial  a  -  .9 

Iteration 

a 

J(a) 

X(l) 

0 

.9 

78.9130445 

-7.7572184 

1 

.8061489 

11.8804011 

-2.7175508 

2 

.7036760 

1.6122215 

-.7432785 

3 

.6046577 

.1689001 

.0158974 

4 

.5327346 

.0090401 

.2877378 

5 

.5034419 

.0000834 

.3602260 

6 

.5000500 

.0000000 

.3677829 

Table  4 . 1 


a  =  1,  r  = 

Iteration 

a 

J(a) 

x(l) 

0 

.2 

.1205106 

.637832 

1 

.8835656 

55.8792983 

-6.443195 

2 

.7876713 

8.3050417 

-2.201167 

3 

. 6845883 

1.0919149 

-.542773 

4 

.5884893 

.1055430 

.090622 

5 

.5236488 

.0045683 

.311016 

6 

.5018820 

.0000246 

.363724 

7 

.5000214 

.0000000 

.367845 

Table  4.2 


Example  4.2.  Here  we  set  a  =  l,  0  =  1,  y  =  4,  and  a  =  .9.  The 
nonhomogeneous  term  is 

,  „10e”^  ^.1 

f(t)  -  -1  -  V(.l)  ^  • 

In  this  case  the  true  solution  is  e 

This  example  contains  a  kernel  that  is  more  singular  than 
that  of  Example  1.  The  results  for  initial  values  of  =  .999 

and  o£q  =  .8  are  given  in  Tables  (4.3)  and  (4.4),  respectively. 

For  comparison, the  true  value  of  x(l)  is  x(l)  =  .3678794.  Note 
again  it  appears  that  a  high  initial  guess  of  a  is  preferable 


57 


to  a  low  one.  Moreover,  for  an  initial  guess  of  =  .75  the 

algorithm  updates  the  parameter  to  a  value  greater  than  1  and 
the  program  stops.  The  convergence  of  the  states  is  shown  in 
Figures  (4.3)  and  (4.4) . 


4,  0=1,  true  a  =  .9,  initial  a  =  .999 

Iteration 

a 

J(a) 

X(l) 

0 

.999 

56.1016515 

7.127846 

1 

.9400160 

4.37204423 

2.223429 

2 

.9075174 

.1077531 

.656523 

3 

.9002876 

.0001462 

.378490 

4 

.9000004 

.0000000 

.367895 

Table  4.3 


a  =  1,  7  = 

4,  |S  =  1,  true 

a 

=  .9,  initial  a 

=  .8 

Iteration 

a 

J(a) 

X(l) 

0 

.8 

7.0279032 

1.895701 

1 

.9676194 

17.3817062 

4.096880 

2 

.9200856 

.8805865 

1.195951 

3 

.9019895 

.0071213 

.441968 

4 

.9000204 

.0000007 

.368632 

5 

.9000000 

.0000000 

.367879 

Table  4 . 4 


Example  4.3.  This  example  has  two  features  that  are  different 
than  the  previous  examples.  In  Section  2  we  assumed  that 
0  >  0,  ensuring  that  the  integral  in  equation  (2.4)  exists. 
This  example  illustrates  that  it  may  be  possible  to  lift  this 
restriction  to  include  0=0,  which  results  in  a  true 
fjagtional  derivative  model.  Also,  for  this  example  x(t)  = 
f*"*  .  Thus  the  true  solution  has  an  unbounded  second 
derivative  at  t  =  0.  Because  integration  methods  based  on 
Simpson's  rule  converge  slowly  for  functions  that  do  not  have 
four  continuous  derivatives,  it  was  necessary  to  increase  N  to 
200  for  this  example  to  gain  accuracy. 

Here  we  set  the  values  a=l,  0  =  0,  r  =  5,  and  a  =  .5. 

The  function  f(t)  is  in  this  case 
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f(t) 


2t 


2.5 


5 


r(.5)y 


The  results  of  the  quasilinearization  algorithm  are  given  in 
Table  (4.5)  for  the  initial  value  of  ot^  =  .9.  The  fact  that  we 

could  only  obtain  a  to  3  correct  digits  is  due  to  the 
inaccuracy  of  the  integration  scheme  used.  Figure  (4.5) 
illustrates  the  convergence  of  the  states  to  the  true 
solution. 


r 


a=l,  r=5,  0=0,  true  a  =  .5,  initial  a  =  .9 


Iteration 

0 

1 

2 

3 

4 

5 

6 


a 

.9 

.7641130 

.6246967 

.5298167 

.5019498 

.5002170 

.5002110 


J(a) 

24.1608906 

3.5153401 

.3481690 

.0127381 

.0000392 

.0000000 

.0000000 


X(l) 

5.563632 

2.711820 

1.531838 

1.101062 

1.005624 

1.000050 

1.000030 


Table  4.5 


59 


References 


1.  R.  L.  Bagley  and  P.  J.  Torvik,  Fractional  calculus  —  a 
different  approach  to  the  analysis  of  viscoelastically  damped 
structures,  AIAA  Journal  21  (1983),  741-748. 

2.  H.  T.  Banks,  R.  H.  Fabiano,  and  Y.  Wang,  Estimation  of 
Boltzmann  damping  coefficients  in  beam  models,  LCDS/CCS  Report 
No.  88-13,  Brown  University,  Providence,  1988. 

3.  H.  T.  Banks  and  G.  M.  Groome,  Jr.,  Convergence  theorems 
for  parameter  estimation  by  quasi linearization,  J.  Math.  Anal. 
Appl.  42  (1973),  91-109. 

4.  H.  T.  Banks  and  K.  Kunisch,  Estimation  Techniques  for 
Distributed  Parameter  Systems,  Birkhauser,  1989. 

5.  V.  Barbu  and  S.  I.  Grossman,  Asymptotic  behaviour  of 
linear  integrodifferential  systems,  Trans.  Amer.  Math.  Soc. 

173  (1972) ,  277-288. 

6.  D.  W.  Brewer,  The  differentiability  ivith  respect  to  a 
parameter  of  the  solution  of  a  linear  abstract  Cauchy  problem, 
SIAM  J.  Math.  Anal.  13  (1982),  607-620. 

7.  D.  W.  Brewer,  A  parameter  dependence  problem  in  functional 
differential  equations.  Lecture  Notes  in  Pure  and  Applied 
Mathematics,  81  (1982),  Marcel  Dekker,  Inc.,  187-195. 

8.  D.  W.  Brewer,  J.  A.  Burns,  and  E.  M.  Cliff,  Parameter 
identification  for  an  abstract  Cauchy  problem  by 

quasi linearization,  ICASE  Report  No.  89-75,  1989. 

9.  J.  A.  Burns  and  R.  H.  Fabiano,  Modelling  and  approximation 
for  a  viscoelastic  control  problem,  Proc.  Third  International 
Conference  on  Control  and  Identification  of  Distributed 
Systems,  Vorau,  July  1986,  23-39. 

10.  J.  A.  Burns  and  T.  L.  Herdman,  Adjoint  semigroup  theory 
for  a  Volterra  integrodifferential  system.  Bull.  Amer.  Math. 
Soc.  81  (1975),  1099-1102. 

11.  J.  A.  Burns  and  T.  L.  Herdman,  Adjoint  semigroup  theory 
for  a  class  of  functional  differential  equations,  SIAM  J. 

Math.  Anal.  7  (1976),  729-745. 

12.  W.  Desch  and  R.  K.  Miller,  Exponential  stabilization  of 
Volterra  integral  equations  with  singular  kernels,  J.  Integral 
Equations  and  Applications  1  (1988),  397-433. 


60 


13.  Y.  C.  Fung,  Foundations  of  Solid  Mechanics,  Prentice-Hall, 
Inc.,  1965. 

14.  K.  B.  Hannsgen  and  R.  L.  Wheeler,  Time  delays  and  boundary 
feedback  stabilization  in  one-dimensional  viscoelasticity , 
Proc.  Third  International  Conference  on  Control  and 
Identification  of  Distributed  Parameter  Systems,  Vorau,  July 
1986,  136-152. 

15.  W.  J.  Hrusa,  J.  A.  Nohel,  and  M.  Renardy,  Initial  value 
problems  in  viscoelasticity,  Appl.  Mech.  Rev.  41  (1988), 
371-378. 

16.  P.  Linz,  Analytical  and  Numerical  Methods  for  Volterra 
Equations,  SIAM  Studies  in  Applied  Mathematics,  1985. 

17.  M.  Renardy,  W.  J.  Hrusa,  and  J.  A.  Nohel,  Mathematical 
Problems  in  Viscoelasticity,  Pitman  Monographs  and  Surveys  in 
Pure  and  Applied  Mathematics  35,  Longman  Scientific  and 
Technical,  Burnt  Mill,  1987. 

18.  P.  J.  Torvik  and  R.  L.  Bagley,  On  the  appearance  of  the 
fractional  derivative  in  the  behaviour  of  real  materials ,  J. 
Appl.  Mech.  51  (1984),  294-298. 


61 


axis 


lire  4.3 


axis 


A  DIRECT  METHOD  FOR  PARAMETER 
ESTIMATION  IN  DISTRIBUTED  SYSTEMS^ 
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and 
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Fayetteville,  Arkansas  72701 


Abstract.  We  consider  identification  of  parameters  in  a 
partial  differential  equation  modeling  the  longitudinal 
vibrations  of  a  viscoelastic  beam.  A  semi-discrete 
approximation  of  this  model  gir es  rise  to  a  Volterra 
inuegrc-differential  system  with  a  weakly  singular  kernel. 
Such  kei.nels  arise  in  fractional  derivative  damping  models  of 
viscoelastic  matei-ials.  A  quasilinearization  method  is  used 
to  identify  damping  parameters  in  the  system.  We  present  the 
results  of  numerical  experiments  using  a  Galerkin 
approximation  in  space  and  an  approximation  in  time  which  is 
adapted  to  weakly  singular  systems  whose  kernels  are  not 
continuous  at  the  origin. 


1.  Introduction.  In  this  paper  we  consider  numerical  methods 
for  the  identification  of  parameters  in  an  idealized  physical 
model  for  the  longitudal  motions  of  a  uniform  bar  fixed  at 
both  ends  with  Boltzmann  type  damping.  The  governing  equation 
is  ([10],  [15]) 


(1.1) 


^  f  9  ) 

pu^^(x,t)  =  —  |Eu^(x,t)  +  ^  J  g(t-s)u^(x,s)ds  J- 


+  f(x,t),  0  <  X  <  1,  t  >  0, 


with  boundary  conditions  u(0,t)  =  0,  u(l,t)  =  0, 


^This  research  was  supported  by  the  Air  Force  Office  of 
Scientific  Research  under  grant  AFOSR-89-0472 . 


67 


and  initial  conditions 


u(x,0)  =  d(x) 


u^.(x,0)  =  v(x)  . 


Here,  u(x,t)  represents  the  axial  displacement  of  position  x 
at  time  t,  p  is  the  density  of  the  material,  E  a  stiffness 
parameter,  f(x,t)  a  forcing  function,  and 


g(s,q) 


s  >  0, 


represents  a  fractional  derivative  damping  term  modified  to 
have  exponential  decay  [ip.  Here  r(-)  denotes  the  gamma 
function  and  q  =  (o£,0)  «  with  0  s  a  <  1  and  0  >  0.  Such 
kernels  ari  e  in  the  study  of  fractional  derivative  models  of 
viscoelastic  structures.  We  refer  the  reader  to  [13],  [19], 
[14],  [16],  and  in  particular  to  [18]  and  the  extensive 
bibliography  therein  for  a  discussion  of  viscoelastic  models 
as  they  relate  to  weakly  singular  kernels. 

In  this  paper  we  consider  the  identification  of  a  and  0. 
Much  of  the  groundwork  for  this  paper  appears  in  [9]  and  [8] 
where  the  identification  of  scalar  systems  is  considered  and 
in  [6]  and  [7]  where  the  theoretical  framework  used  here  is 
established.  Torvik  and  Bagley  [1],  [19]  have  estimated  the 
parameter  a  in  the  Laplace  transform  domain.  Banks,  et.  al. , 
[2],  [4]  have  identified  parameters  corresponding  to  $  and  r 
in  a  similar  but  different  model  assuming  a  is  known. 

A  Galerkin  approximation  is  used  to  approximate  (1.1)  by 
a  Volterra  equation  with  weakly  singular  kernel.  In  the 
example  presented  below,  this  semi-discrete  model  is  solved 
using  a  product  integration  method  based  on  Simpson's  rule. 
Product  integration  methods  using  polynomial  collocation  are 
well  suited  for  equations  of  the  form  (1.1)  if  the  solution  is 
analytic  on  the  entire  interval  of  interest.  However,  for  the 
kernel  of  interest  in  this  paper,  the  solution  is  not  analytic 
at  t  =  0  if  f(t)  is  smooth  [17]. 

If  one  applies  a  Galerkin  scheme  to  equation  (1.2),  one 
obtains  a  system  of  second-order  integro-dif ferential 
equations.  Integrating  a  scalar  form  of  this  system  yields  an 
equation  of  the  form  (See  [9]  for  details.) 

a  r  x(s)ds  +  [  g(t-s)x(s)ds  +  f(t), 

J  0  -^0 


A  standard  assumption  in  viscoelasticity  [14]  is  that  the 
material  is  in  an  unstrained  state  for  time  t  <  0.  This  would 
correspo^id  to  u(x,s)  =  0  for  s  <  0  in  equation  (1.1).  Define 

z(t)  =  I  x(s)ds.  It  follows  then  that  x(s)  =  0  and  z(s)  =  0 
•'O 
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for  s  <  0.  If  we  define  w(t)  =  col (x(t)  ,  z  (t) ) ,  then  w(t) 
satisfies 


(1.3) 


(  w(t)  =  Mw(t)  +J  K{t-s)w(s)ds  +  F(t), 


f  X  1 

'  0  ’ 

w(0)  = 

U 

,  w(s)  = 

,  S  <  0, 

0 

.  0  , 

(?o). 

,  K(S)  =  ( 

o  o 

^ _ > 

,  and  F(t)  =  [ 

2.  The  algorithm.  In  this  section  we  state  for  reference  the 
quasilinearization  algorithm  which  is  defined  in  an  abstract 
semigroup  formulation  of  equation  (1.3).  Details  of  this 
construction  and  related  smoothness  results  may  be  found  in 
[7]  and  [9].  We  consider  equation  (1.3)  in  the  form 


(2.1) 


w(t)  =  Mw(t)  +  J  K(--s,q)w(t+s)ds  +  F(t),  t  >  0, 

"^CO 

^  W(0)  =  7),  W(S)  =  <P(S),  S  <  0, 


where  t) 


€  M  = 


a 

0 


and 


K(s,q) 


ye  f  1  ® 

r(l-a)s“  ^  °  ° 


s  >  0. 


We  are  interested  in  dependence  on  the  parameter  q  = 

(a,P)  where  ^  >  0  and  0  s  a  <  i.  By  a  solution  of  (2.1)  we 
mean  a  function  w:  (-w,™)  -♦  such  that  w  is  absolutely 
continuous  (A.C.)  on  [0,<»)  and  satisfies  the  integral  equation 
a.e.  on  [0,®),  w(0)  =  tj,  and  w(s)  =  p(s)  a.e.  on  (-w,0]. 

Our  semigroup  formulation  follows  the  construction  in  [5] 
as  further  de^voloped  in  [11]  and  [12].  Define  the  product 
space  X  =  W  X  L^(-«,0)  with  norm 

ll(T?,^)ll_  =  |77l  +  HiPtLi,  nx  •  Then  it  can  be  shown  that  the 

A  JLi  \  »  u  } 

solution  of  (2.1)  is  given  by  the  variation  of  constants 
formula 


(w(t),w^)  =  S(t,q)  (7],«>)  +  Q(t,q) 

where 
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(2.2) 


Q(t,q)  =  S(t-s,q) (F(s) ,  0)ds. 

•^0 

and  S(t,q)  is  a  CQ-semiqroup  on  X  corresponding  to  the 

homogeneous  version  of  (2.1). 

We  are  now  in  a  position  to  define  a  parameter  estimation 
algorithm  based  on  quasilinearization  and  state  some  local 
convergence  results.  For  later  adaptation  we  develop  the 
algorithm  for  the  case  q  e  P  c  r"  with  canonical  basis 

i  =  1,2,.  .  .,n.  In  Section  3  the  algorithm  is  applied  in 

the  cases  n  =  l  and  n  =  2.  The  definitions  and  theorems 
stated  here  may  also  be  found  in  [7],  but  are  included  here 
for  immediate  reference. 

Let  Yq  =  (Tlj^p)  €  X  and  q  e  P.  Let  C  be  a  bounded  linear 
mapping  from  X  into  a  finite-dimensional  space  R^,  and  define 

y(t,q)  =  CtS(t,q)yQ  +  Q(t,q)3. 

The  parameter  estimation  algorithm  is  related  to  the  following 
optimization  problem. 

I 

Problem  2.1,  Let  YjCR,  3  =1,2,.  .  .,mbe  data  values 
taken  at  times  tj6[0,T],  j  =  l,2,.  .  .,m,  respectively.  For 
q  €  P  define  the  quadratic  cost  function 

m 

J(q)  =  I  ly(tj,q)  - 

j=l 

Find  q*  €  P  such  that  J(q*)  s  J(q)  for  all  q  €  P. 

The  quasilinearization  algorithm  method  defines  a 
recursive  algorithm  whose  fixed  point  is  a  local  solution  of 
Problem  2.1.  A  more  complete  exposition  is  given  in  [3]. 

Given  an  initial  guess  q^  e  P  define 

^k+1  ^  '  k  =  0,1,2,.  .  . 

where 

^(q)  =  ^  -  [D(q)]"'b(q) 

D(q)  =  'I  M^t.,q)M(tj,q) 
m 

b(q)  =  M^(t^,q)  [y(tj,q)  -  y  ^  ] 

j=l  i 

and  the  matrix  M(t,q)  has  its  i  column  M  (t,q)  given  by 
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M^(t,q)  =  CDq[S(t,q)yQ  +  Q(t,q)]e^,  i  =  1,2,.  .  .n. 

The  following  theorems  are  typical  of  quasilinearization 
methods.  Their  proofs  may  be  found  in  [7].  We  obtain 
super linear  convergence  when  there  is  an  exact  fit  to  data 
(Theorem  2.1)  and  linear  convergence  in  the  presence  of  error 
(Theorem  2.2). 

Theorem  2.1.  Suppose  the  conditions  of  the  previous  section 
are  satisfied.  Moreover,  assume  [D(q)]”  exists,  ((q*)  =  q*, 
and  J(q*)  =  0.  Then  for  every  c  >  0  there  exists  S  >  0  such 
that 


U(q)  -  ((q*) 1  s  clq  -  q*| 

for  Iq  -  q*|  ^  5.  In  particular,  there  is  a  neighborhood  ll  of 
q*  such  that  qj^  q*  as  k  »  whenever  q^  e  It. 

The  following  theorem  does  not  regaire  an  exact  fit  to 
data  but  does  place  some  technical  restrictions  on  the 
behaviour  of  the  matrix  M(t,q)  near  q*.  Note  that  under  the 
conditions  of  Theorem  2.2  there  exists  number  6  >  0  such  that 
for  0  <  5  <  5  there  exists  a  constant  K(6)  such  that 

m 

I  IM^tj,q)  -  M^(tj,q*)|  s  K(S)  Iq  -  q*| 

j=l 

for  Iq  -  q*|  s  5.  Let  K*  =  lim  sup  K(5)  and  define 

6  0 

\*  =  K*|D(q*)~*l  max|w(t.,q*)  -  w .  1 . 

j  ^  ^ 

Theorem  2.2.  Suppose  the  conditions  of  t^e  previous  section 
are  satisfied.  Moreover,  assume  tD(q*)]"*  exists  and 
^(q*)  =  q**  Let  X*  be  defined  as  above  and  assume  X*  <  1. 

Then  there  exists  5*  >  0  such  that 

l^(q)  “  #(q*) I  ^  A*lq  -  q*l 

for  Iq  -  q*|  £  5*.  In  particular,  qj^  q*  as  k  -»  «  whenever 
Iqp  -  q*l  =  S*. 


3.  Numerical  results.  In  this  section  we  present  numerical 
results  for  parameter  estimation  in  a  semi-discrete 
approximation  of  the  partial  differential  equation 
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I 


t  -/3(t-s) 

- oT  fl(x,t) 


r(l-a)  at  •'0  (t-s) 

with  boundary  conditions  u(0,t)  =  0,  u(l,t)  =  0,  and 
initial  conditions  u(x,0)  =  Uq(x),  u^(x,0)  =  Uj^(x). 

Integrating  with  respect  to  t  one  obtains 

j"  pu^(x,t)  =  EJ  Uj^^(x,s)ds  + 

t  -0(t-s) 

- - -  - ^  Uyy(x,s)ds  +  f,(x,t) 

r(l-a)  Jq  (t-s)“  ^ 


where  f2(x,t)  =  |  f^(x,s)ds  +  pUj^(x) 


We  apply  a  Galerkin  approximation  in  which  the  interval  [0,1] 
is  divided  into  N  equal  parts  and  the  homogeneous  boundary 
conditions  allow  us  to  approximate  the  solution  by  a  function 

N 

of  the  form  u(x,t)  =  V  a.(t)^.(x)  where  <p.  is  a  cubic  spline 

jko  3  J  J 

basis  element.  Substituting  in  the  equation,  taking  an  inner 
product  with  a  basis  element  and  integrating  by  parts 

yields  a  Volterra  equation  in  t  of  the  form 

I  pA*^v(t)  =  -ED**J  v(s)ds  - 

,t  -0(t-s) 


(3.1) 


r(l-a) 


0  (t-s) 


—  v(s)ds  +  f^Ct) 


N  N 

where  A  and  D  are  N+1  x  N+1  matrices  depending  on  inner 
products  of  and  (p'^. 

The  quasilinearization  algorithm  requires  that  we  solve 
(3.1)  along  with  its  sensitivity  equations  which  have  the  form 
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=■  nW  Jo 


(3.2) 


t  ^-^(t-s) 


(t-s) 


a  a 


X  (s)ds 


r(l-a)^  ■’o 


(1-a) 


(t-s) 


a 


,  y  ln(t-s)e~^^^'‘^^  xfslds 

Jn  (t-s)“ 


=  0, 


f  x^(t) 


(3.3) 


■  *  rrW  I 


r( 


^Jo 


0  (t-s) 

t  -0(t-S) 


t  -0(t-s) 


(t-s) 


x^(0) 


=  0, 


(t-s) 


a 


x(s) ds 


where  ^  and  =  -r^.  The  zero  initial  condition 

OC  a<X  p  dp 

reflects  the  fact  that  x(0)  =  is  independent  of  a  and  0. 

It  is  important  to  note  that  (3.2)  and  (3.3)  also  have  weakly 
singular  kernels. 

The  impleraentaton  of  the  identification  scheme  begins 
with  an  initial  guess  for  q=  (a,/3).  Equations  (3.1)-(3.3)  are 
integrated  using  this  initial  value,  then  x(t),  x^(t),and 

x^(t)  are  used  to  give  an  updated  estimate  of  the  parameter 

using  the  quasilinearization  algorithm. 

As  a  specific  example  consider  the  partial  differential 
equation 


f  pu^^(x,t)  =  Eu^^(x,t) 


J 


t^-0(t-s) 


r(i-a)  at 

with  boundary  conditions 


a 


0  (t-s) 

u(0,t)  =  0,  u(l,t) 
t 


Uxx(XfS)ds  +  f(x,t) 


0,  and  initial 


conditions  u(x,0)  =  0,  u.  (x,0)  =  47rsin(2rrx)  ,  where  p  =  l, 


E  =  l,  y  =  l,  a  =  1/2,  and  0  =  0.  The  exact  solution  is 
u(x,t)  =  sin(27rx)sin(47it)  with  the  forcing  function  chosen  as 
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The  numerical  scheme  uses  N  =  8  and  h  =  .02.  The  parameter 
estimator  uses  70  data  points  taken  from  the  exact  form  of 
u^(x,t) .  The  results  of  an  application  of  the 

quasilinearization  algorithm  to  the  estimation  of  the 
parameter  a  are  given  in  the  table  below. 


Iteration  a  J(a) 

0  .9  1703.2400606 

1  0.0783169  10811.4763257 

2  0.2310107  4486.6256214 

3  0.3658041  1012.4621853 

4  0.4729954  32.7112404 

5  0.4990789  0.0523379 

6  0.5002075  0.0001728 

7  0.5002095  0.0001726 

8  0.5002095  0.0001726 


The  table  below  shows  the  results  of  the  algorithm  for 
identifying  both  the  parameter  a  and  |3. 


Iteration 

0 

1 

2 

3 

4 

5 

6 

7 

8 


a 

0.2000000 

0.2986892 

0.5281903 

0.4674188 

0.5135026 

0.5203233 

0,5029341 

0.5001956 

0.5002027 


0.5000000 

6.5002076 

7.7538270 

-4.2094682 

-3.4408819 

-2.3067455 

-0.1314789 

0.0032037 

0.0009508 


5076.0182126 

1394.6849127 

273.4823741 

469.6983381 

70.7243517 

26.1540142 

0.2291736 

0.0001997 

0.0001674 
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